import pandas as pdAlgoritmo para la estimación de la turbidez sobre el Río Paraná
Este sitio web contiene información sobre la estimación de la turbidez por teledetección en la cuenca media del Río Paraná. La turbidez es uno parámetros de interés dentro proyecto Estimar indicadores de calidad de agua en la cuenca media del río Paraná para el desarrollo de un algoritmo mediante técnicas de teledetección satelital (MSECRE0008604), desarrollado por el Grupo de Investigación Sobre Temas Ambientales y Químicos (GISTAQ) de la Universidad Tecnológica Nacional Facultad Regional Resistencia (UTN-FRRe).
Se utilizarán imágenes del satélite Sentinel-2 con corrección automática, de las cuales se obtiene la reflectancia de superficie del agua. Se buscará la relación entre la reflectancia y la turbidez por métodos de regresión tradionales y machine learning. Una vez obtenido el algoritmo que relacione ambas propiedades, se desarrollaran mapas de distribución espacial.
GISTAQ, UTN, FRRe, Algoritmo, Turbidez, Machine learning, Teledetección
1 Turbidez
La turbidez se refiere a la opacidad o falta de claridad en un líquido provocada por la presencia de partículas suspendidas. Este fenómeno es un indicador clave en el monitoreo de la calidad del agua y su influencia en diferentes ecosistemas es significativa.
La turbidez es un indicador de la calidad del agua, reflejando la presencia de partículas en suspensión. Su medición es crucial para garantizar la potabilidad del agua y la salud de los ecosistemas acuáticos. Este fenómeno puede ser resultado de diversas causas, como la erosión del suelo, la actividad biológica y la contaminación. La comprensión de la turbidez y su impacto es esencial para la gestión de recursos hídricos y la protección del medio ambiente.
La turbidez viene determinada por la dispersión de la luz causada por la materia suspendida en el agua, se obtiene normalmente mediante un turbidímetro, que proporciona medidas en Nephelometric Turbidity Unit (NTU) y mide la dispersión de un rayo de luz en el agua a 90º de la luz incidente [1].
Muchas propiedades, como la clorofila-a (Chl-a), sólidos suspendidos totales (SST) y la materia orgánica disuelta coloreada (CDOM), se utilizan a menudo como indicadores del estado del agua. Estos constituyentes del agua a su vez son responsables de la turbidez.
Existe una fuerte correlación entre turbidez y sólidos suspendidos totales, por lo que se puede estimar SST a partir de la turbidez. Por lo general, es una relación directa, a mayor concentración de SST mayor turbidez.
Existe una relación inversa entre la Turbidez y la profundidad del disco de Secchi (a valores bajos de secchi mayor turbidez), por lo que también se puede estimar turbidez a partir de mediciones de disco de secchi.
1.1 Métodos tradicionales
| Ecuación | Bandas (nm) | Métricas | Aguas | Plataforma | Referencia |
|---|---|---|---|---|---|
| \(1.559e^{35.533 \cdot B03} \\ 1.879e^{37.745(B03 \cdot B5)/(B04+B12)}\) | B03, B04, B05, B12 | \(R^{2}\), RMSE, MAE | Lago1 | Sentinel-2 | [2] |
| \(2677.2 \cdot B04^{1.856}\) | B04 | \(R^{2}\), RMSE, Bias | Interiores variadas2 | Landsat-8 | [3] |
| \(969-1.5468 \cdot R_{1200nm}+2.07 \frac{B8A}{B02}\) | B02, B8A, 1200nm | IOA, SI, RMSE, MAE | Río3 | Landsat-8 | [4] |
| \(y=-1.1+5.8 \frac{B02}{B04} \\ y=3.896-4.186 \frac{B02}{B03}\) | B02, B03, B04 | \(R^{2}\), RMSE | Río4 | Landsat-8 | [5] |
| \(y=37661 \cdot B8A^{2}+1845 \cdot B8A <br> y=531.5- \frac{B04}{0.88}\) | B04, B8A | \(R^{2}\), RMSE, MAPE | Estuario5 | Pléiades | [6] |
Múltiples modelos (lineal, logaritmos, inversa, cuadrática, exponencial, potencial) y plataformas (Sentinel-2, Landsat-5 y Landsat-8) emplean el cociente de bandas B04/B03 [7].
Modelos de estimación a partir de Sentinel-2 y Landsat-8 utilizan regresiones lineales, cuadráticas y logarítmicas empleando B02, B03, B04, B01 (con menos apariciones) y cocientes entre éstas [8].
1.2 Métodos de aprendizaje automático
El aprendizaje automático es un subconjunto de la inteligencia artificial que permite que un sistema aprenda y mejore de forma autónoma, sin necesidad de una programación explícita, a través del análisis de grandes cantidades de datos. El aprendizaje automático permite que los sistemas informáticos se ajusten y mejoren continuamente a medida que acumulan más “experiencias”. Por lo tanto, el rendimiento de estos sistemas puede mejorar si se proporcionan conjuntos de datos más grandes y variados para su procesamiento.
Cuando se entrenan modelos de machine learning, cada conjunto de datos y cada modelo necesitan un conjunto diferente de “hiperparámetros”. Los hiperparámetros son variables de configuración externa que se utilizan para administrar el entrenamiento de modelos de machine learning. Controlan de forma directa la estructura, funciones y rendimiento de los modelos. Los hiperparámetros son los parámetros de un modelo de aprendizaje automático, que no se aprenden durante el entrenamiento, sino que se establecen antes de que comience.
El “ajuste de hiperparámetros” permite modificar el rendimiento del modelo para lograr resultados óptimos. Este proceso es una parte fundamental del machine learning. El ajuste de hiperparámetros puede ser manual o automático. A pesar de que el ajuste manual es lento y tedioso, permite entender mejor cómo afectan al modelo las ponderaciones de los hiperparámetros. El proceso de ajuste de hiperparámetros es iterativo, y debe probar diferentes combinaciones de parámetros y valores.
En el aprendizaje automático es importante utilizar técnicas de “validación cruzada” , de modo que el modelo no se centre únicamente en una única porción de sus datos. La validación cruzada o cross-validation es una técnica utilizada para evaluar los resultados de un análisis estadístico y garantizar que son independientes de la partición entre datos de entrenamiento y prueba. La idea básica de la validación cruzada es dividir los datos en conjuntos de entrenamiento y validación, y luego entrenar el modelo en el conjunto de entrenamiento y evaluar su rendimiento en el conjunto de validación. Este proceso se repite varias veces, con diferentes subconjuntos de los datos utilizados para el entrenamiento y la validación, y se calcula el rendimiento promedio.
En los procesos de machine learning supervisado se utilizan diversos algoritmos y técnicas de cálculo, generalmente calculados mediante el uso de programas como R o Python.
Dependiendo del tipo de datos que se usen para el entrenamiento, será de modelo de aprendizaje automático que se use. A grandes rasgos, existen tres tipos de modelos que se usan en el aprendizaje automático: aprendizaje supervisado , no supervisado y por refuerzo.
Consultando el trabajo de otros investigadores, se observa que utilizan principalmente el aprendizaje automático supervisado. Este tipo aprendizaje supervisado utiliza un conjunto de entrenamiento para enseñar a los modelos a producir el resultado deseado. Este conjunto de datos de entrenamiento incluye entradas y salidas correctas, que permiten al modelo aprender con el tiempo. El algoritmo mide su precisión a través de la función de pérdida, ajustando hasta que el error se haya minimizado lo suficiente.
Yang Zhe y otros, utilizaron como datos de entrada la reflectancia de superficie y datos de salida la turbidez, utilizaron los modelos SVR (support vector regression), random forest (RF) y eXtreme Gradiente Boostring (XGBoost). Los hiperparámetros de cada modelo se determinaron mediante una búsqueda en cuadrícula de validación cruzada en Scikit-Learn de Python [9].
Ma Yue y otros, utilizaron varios modelos de aprendizaje automático, usaron Python 3.7 tanto para la predicción de la turbidez del agua y la optimización de la los hiperparámetros [2].
Zhao y otros probaron 14 modelos de machine learning en un estanque de peces con un dispositivo de construction propia, de los cuales ETR, Bagging, RFR, and ABR son los que presentaron un mejor desempeño en la estimación de la turbidez. Los algoritmos se implementaron utilizando Python 3.6 y bibliotecas de aprendizaje scikit [10].
| Modelo de machine learning | Cuerpo de agua | Métricas | Plataforma | Referencia |
|---|---|---|---|---|
| SVR, ELM ,BP ,CART ,GBT ,RF ,KNN | Lagos | RMSE, \(R^{2}\), MAE | Sentinel-MSI | [2] |
| eXtreme Gradient Boosting (XGBoost), support vector regression (SVR), random forest (RF) | Lago | RMSE, \(R^{2}\), MAPE | Sentinel-2A/B y Landsat-8/9 | [9] |
| linear regression (LR), ridge regression (RR), least absolute shrinkage and selection operator regression(LASSO), elastic net regression (ENR), k-nearest neighbor regression (KNN), Gaussian process regression (GPR), decision tree regression (DTR), support vector regression (SVR), multilayer perceptron regression (MLP), adaptive boosting regression (ABR), gradient boosting regression (GBR), bootstrap aggregating regression (Bagging), random forest regression (RFR), and extreme tree regression (ETR) | Estanque de peces | MAE, MRSE, MAPE, \(R^{2}\), RE, Acc | Propia | [10] |
2 Procesamiento de datos
Para el procesamiento de los datos se utilizará la librería pandas de Python.
En el proyecto tenemos dos archivos .csv que contienen los datos:
-base_de_datos_lab.csv → contiene resultados de laboratorio
-base_de_datos_gis.csv → contiene datos espectrales
Importamos la librería pandas para usarla, la nombramos como “pd” para simplificar
Leemos los archivos .csv por separado y definimos dos DataFrame
Un DataFrame es basicamente una tabla, donde la información se organiza en filas y columnas. Los datos de la misma columna contienen el mismo tipo de datos, pandas agrega por defecto un “índice” que nos ayuda a identificar una fila en particular.
Con la función pd.read_csv le idicamos a pandas que queremos leer archivos .csv.
df1_lab = pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_lab.csv")
df2_gis = pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_gis_acolite.csv")df1_lab DataFrame de datos provenientes del laboratorio.
df2_gis DataFrame de datos espectrales provenientes del sensor MSI de Sentinel-2.
Nota: Se debió colocar la “r” delante de la dirección para que lea los archivos.
Video de YouTube ¿Qué es un DataFrame?
Verificamos que los datos se han leído y se crearon correctamente ambos DataFrame por separado, con print
df1_lab.head()
df2_gis.head()| punto | banda | reflect | fecha | longitud | latitud | |
|---|---|---|---|---|---|---|
| 0 | 1 | B01 | 0.041367 | 2023-05-11 | -58.868047 | -27.464687 |
| 1 | 1 | B02 | 0.046669 | 2023-05-11 | -58.868047 | -27.464687 |
| 2 | 1 | B03 | 0.085630 | 2023-05-11 | -58.868047 | -27.464687 |
| 3 | 1 | B04 | 0.128799 | 2023-05-11 | -58.868047 | -27.464687 |
| 4 | 1 | B05 | 0.128704 | 2023-05-11 | -58.868047 | -27.464687 |
Nota: La primer columna (donde se ven los valores 0,1,2,3,4, …) es el índice que agrega pandas por defecto al DataFrame, esa columna no forma parte del csv original.
Combinamos ambos DataFrame para tener los datos en una única base de datos.
Lo defimos como df_combinado, para esto utilizamos la función pd.merge para realizar la combinación.
df_combinado = pd.merge(df1_lab, df2_gis, on=['latitud', 'longitud','fecha'], how='inner')on=[‘latitud’, ‘longitud’,‘fecha’] Especifica las columnas por las cuales se unirán los dos DataFrames. En este caso, por coincidencias exactas en latitud, longitud y fecha.
how=‘inner’ Es el tipo de unión, significa que solo se conservarán las filas que tengan coincidencias en ambos DataFrames en las columnas mencionadas.
Verificamos que la combinacion se haya realizado correctamente
df_combinado.head()| fecha | longitud | latitud | param | valor | punto | banda | reflect | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2023-05-11 | -58.868047 | -27.464687 | ph | 6.8 | 1 | B01 | 0.041367 |
| 1 | 2023-05-11 | -58.868047 | -27.464687 | ph | 6.8 | 1 | B02 | 0.046669 |
| 2 | 2023-05-11 | -58.868047 | -27.464687 | ph | 6.8 | 1 | B03 | 0.085630 |
| 3 | 2023-05-11 | -58.868047 | -27.464687 | ph | 6.8 | 1 | B04 | 0.128799 |
| 4 | 2023-05-11 | -58.868047 | -27.464687 | ph | 6.8 | 1 | B05 | 0.128704 |
Filtramos la turbidez del DataFrame
Del DataFrame combinado, solo nos interesa las filas que contengan los valores de turbidez, ya que es la propiedad de estudio en este sitio web. Por ello nos quedamos con las filas en donde la columna param sea igual a turb.
Creamos un nuevo DataFrame a partir de df_combinado.
df_turbidez = df_combinado[(df_combinado['param'] == 'turb')]Verificamos
df_turbidez.head()| fecha | longitud | latitud | param | valor | punto | banda | reflect | |
|---|---|---|---|---|---|---|---|---|
| 33 | 2023-05-11 | -58.868047 | -27.464687 | turb | 185.5 | 1 | B01 | 0.041367 |
| 34 | 2023-05-11 | -58.868047 | -27.464687 | turb | 185.5 | 1 | B02 | 0.046669 |
| 35 | 2023-05-11 | -58.868047 | -27.464687 | turb | 185.5 | 1 | B03 | 0.085630 |
| 36 | 2023-05-11 | -58.868047 | -27.464687 | turb | 185.5 | 1 | B04 | 0.128799 |
| 37 | 2023-05-11 | -58.868047 | -27.464687 | turb | 185.5 | 1 | B05 | 0.128704 |
Nota: Se han eliminado las filas de ph, cond, sol_sus, etc. Se conserva el índice del DataFrame original.
Eliminamos las columnas que no nos interesan
Creamos un nuevo DataFrame a partir de df_turbidez. Con la función .drop especificamos que columnas queremos eliminar.
df_turbidez_banda = df_turbidez.drop(columns=['longitud','latitud','punto','fecha','param'])
df_turbidez_banda.head()| valor | banda | reflect | |
|---|---|---|---|
| 33 | 185.5 | B01 | 0.041367 |
| 34 | 185.5 | B02 | 0.046669 |
| 35 | 185.5 | B03 | 0.085630 |
| 36 | 185.5 | B04 | 0.128799 |
| 37 | 185.5 | B05 | 0.128704 |
Nota: Solo nos quedamos con las columnas valor,banda y reflect.
Cambiamos el nombre de la columna “valor” por el de “turbidez” Como los valores de esa columna son de turbidez, directamente cambiamos el nombre de la columna con la función .rename
df_turbidez_banda.rename(columns={'valor': 'turbidez'}, inplace=True)
df_turbidez_banda.head()| turbidez | banda | reflect | |
|---|---|---|---|
| 33 | 185.5 | B01 | 0.041367 |
| 34 | 185.5 | B02 | 0.046669 |
| 35 | 185.5 | B03 | 0.085630 |
| 36 | 185.5 | B04 | 0.128799 |
| 37 | 185.5 | B05 | 0.128704 |
Creamos la tabla final
Usamos la función .pivot_table para que las columnas sean la turbidez y las distintas bandas de Sentinel-2 (B01, B02 , B03…)
df_final = df_turbidez_banda.pivot_table(
index='turbidez',
columns='banda',
values='reflect',
)¿Qué significa cada término en el argumento de la función pivot_table?
index=‘turbidez’ → Las filas serán los valores únicos de ‘turbidez’
columns=‘banda’ → Las columnas serán los valores únicos de ‘banda’ (B01, B02…)
values=‘reflect’ → El contenido de la tabla será lo que haya en la columna ‘reflect’
Creamos archivo final con los valores de turbidez y reflectancia de cada banda
-Forma 1 Creamos un .csv con los datos que nos interesan con la función .to_csv
df_final.to_csv('Turbidez_FINAL.csv', index=True)IMPORTANTE: index=True , debe ser así para que el índice que definimos al pivotar sea una columna visible en el archivo csv.
Recordemos que la columna index es solo visible solo por la librería pandas, Si guardamos con index=False, se omite y no se guarda en el csv.
-Forma 2 Si decidimos poner index=False tenemos que usar una función adicional antes de exportar, debido a que la turbidez está en el índice, no como una columna.
Por lo tanto, luego de hacer el pivot se debe agregar una línea de código.
df_final = df_turbidez_banda.pivot_table(
index='turbidez',
columns='banda',
values='reflect',
)
df_final = df_final.reset_index() #Línea de código que se debe agregarreset_index() convierte el índice definido previamente como “turbidez” en una columna normal y reemplaza el índice por uno numérico estándar (0, 1, 2…). que es el predeterminado por pandas.
Finalmente exportamos como archivo csv.
df_final.to_csv('Turbidez_FINAL2.csv', index=False)Verificación tabla final
import pandas as pd
df = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
df.head()| turbidez | B01 | B02 | B03 | B04 | B05 | B06 | B07 | B08 | B11 | B12 | B8A | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2.0 | 0.030264 | 0.035694 | 0.052761 | 0.038619 | 0.028929 | 0.017878 | 0.018247 | 0.014811 | 0.008830 | 0.006162 | 0.015677 |
| 1 | 3.0 | 0.018744 | 0.027450 | 0.044980 | 0.028040 | 0.018781 | 0.007693 | 0.007571 | 0.007623 | 0.001773 | 0.000112 | 0.003305 |
| 2 | 4.0 | 0.029141 | 0.036459 | 0.054631 | 0.039793 | 0.029767 | 0.022395 | 0.022966 | 0.019097 | 0.009430 | 0.007976 | 0.020290 |
| 3 | 4.5 | 0.022306 | 0.027934 | 0.045472 | 0.028727 | 0.020835 | 0.014462 | 0.017479 | 0.013988 | 0.000088 | 0.000361 | 0.013232 |
| 4 | 5.0 | 0.029902 | 0.030043 | 0.047235 | 0.033111 | 0.024060 | 0.011116 | 0.012519 | 0.007006 | 0.001982 | 0.000181 | 0.009074 |
3 Prueba de modelo lineal
Regresión lineal con método de mínimos cuadrados
Reemplamos el código del tutorial de sklearn por nuestros datos.
Debemos importar las funciones necesarias para:
1- Leer el archivo .csv creado previamente , que contiene los datos de turbidez y los valores de reflectancia para cada banda.
2- Para usar el método de mínimos cuadrados
Importamos las funciones
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as pltpandas → para leer los datos
train_test_split → para dividir los datos en entrenamiento y validación
LinearRegression → para crear el modelo de regresión lineal
mean_squared_error y r2_score → Para medir el desempeño del modelo (RMSE y R2)
matplotlib.pyplot → para realizar gráficos
Leemos los datos de interés y los dividimos en entrenamiento y validación.
Lectura de datos, lo hacemos con pd.read_csv
df = pd.read_csv(r"D:\GIT\Turbidez\Turbidez_FINAL.csv")
X = df[['B05']]
y = df['turbidez'] Definimos:
df, un DataFrame que contiene los datos de turbidez y reflectancia;
X para que sea los valores de reflectancia de una banda en particular;
y para que sea la turbidez.
Dividimos en datos para el entrenamiento y validación
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, shuffle=True, random_state=42
)Conjunto de datos para el entrenamiento X_train y X_test
Conjunto de datos para la validación y_train y y_test
Creamos el modelo de regresión
regressor = LinearRegression().fit(X_train, y_train)LinearRegression() Crea el modelo de regresión lineal llamado “regressor”. Este modelo, se entrena con la función .fit a partir de los datos de entrenamiento (X_train, y_train)
Evaluamos el modelo generado a partir de las métricas de desempeño.
y_pred = regressor.predict(X_test)
p_rmse = mean_squared_error(y_test, y_pred)
p_r2 = r2_score(y_test, y_pred)Con la función print visualizamos los valores de las métricas de desempeño
print("RMSE", p_rmse)
print("R2", p_r2)RMSE 1151.59791674541
R2 0.27485564458517464
Para ver la ecuación del modelo de regresión
# Coeficientes (pendientes)
coef = regressor.coef_
# Intercepto (ordenada al origen)
intercept = regressor.intercept_
print(f"La ecuación es: y = {coef[0]:.3f}x + {intercept:.3f}")La ecuación es: y = 1316.107x + -41.516
Visualizamos los resultados comparando el conjunto de entrenamiento y validación.
Código
fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)
#Gráfico de entrenamiento
ax[0].plot(
X_train,
regressor.predict(X_train),
linewidth=3,
color="#17A77E",
label="Modelo",
)
ax[0].scatter(X_train, y_train, label="Entrenamiento", color = "#9D50A6", alpha = .6)
ax[0].set(xlabel="Banda 5", ylabel="Turbidez", title="Conjunto de entrenamiento")
ax[0].legend()
#Gráfico de validación
ax[1].plot(X_test, y_pred, linewidth=3, color="#17A77E", label="Modelo")
ax[1].scatter(X_test, y_test, label="Validación", color = "#9D50A6", alpha = .6)
ax[1].set(xlabel="Banda 5", ylabel="Turbidez", title="Conjunto de validación")
ax[1].legend()
#Ecuación de la recta
coef = regressor.coef_[0]
intercept = regressor.intercept_
equation = f"y = {coef:.2f}x + {intercept:.2f}"
# Mostrar la ecuación en ambos subgráficos (opcionalmente, puedes usar solo uno)
for a in ax:
a.text(0.05, 0.95, equation, transform=a.transAxes,
fontsize=10, verticalalignment='top',
bbox=dict(boxstyle="round", facecolor="white", alpha=0.7))
fig.suptitle("Regresión lineal")
plt.show()4 Pruebas de correlación
Para ser más rigurosos, agregamos más etapas durante el entrenamiento de nuestro modelo lineal.
En esta etapa nos interesa estimar la turbidez por medio de modelos lineales sencillos, para encontrar mejores modelos podemos aplicar transformaciones logarítmicas con el objetivo principal de linealizar relaciones no lineales.
Primeramente se probará un modelo sin aplicar ninguna transformaciones en las variables.
En las siguientes secciones, se aplicará logaritmo en la turbidez, en las bandas espectrales y en en ambas, con el fin de evaluar en rendimiento de los modelos y seleccionar los mejores, según sus métricas de desempeño.
Se utilizará el coeficiente de correlación lineal r, para evaluar la relación lineal entre la turdidez y las distintas bandas (con y sin transformación logarítmica).
Esta medida indica cuanto se relacionan dos variables, puede tomar valores desde -1 a +1:
• +1 correlación lineal perfecta positiva.
• 0 sin correlación.
• -1 correlación lineal perfecta negativa.
4.1 Prueba de correlación C1: (sin aplicar logaritmo)
Importamos pandas para leer los datos
Código
import pandas as pd
Datos= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
#print (Datos.head())Calculamos coeficiente de correlacion lineal “r” entre la turbidez y cada banda con la función .corr de pandas.
Código
bandas = [col for col in Datos.columns if col.startswith('B')]
correlaciones = {}
for banda in bandas:
r = Datos['turbidez'].corr(Datos[banda])
correlaciones[banda] = r
#Creamos un Data Frame.
df_correlaciones1 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones1 = df_correlaciones1.sort_values(by='r', ascending=False).reset_index(drop=True)
df_correlaciones1['Banda'] = "turb vs " + df_correlaciones1['Banda'].astype(str)
df_correlaciones1.rename(columns={'Banda': 'Combinación 1'}, inplace=True)
df_correlaciones1.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv', index=False)
#Creo archivo csv para usarlo para entrenar modelos
Datos.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv", index=False)
df_correlaciones1| Combinación 1 | r | |
|---|---|---|
| 0 | turb vs B05 | 0.887299 |
| 1 | turb vs B06 | 0.876827 |
| 2 | turb vs B08 | 0.869149 |
| 3 | turb vs B07 | 0.867741 |
| 4 | turb vs B04 | 0.827019 |
| 5 | turb vs B8A | 0.788088 |
| 6 | turb vs B03 | 0.673358 |
| 7 | turb vs B02 | 0.614315 |
| 8 | turb vs B01 | 0.611305 |
| 9 | turb vs B12 | 0.429200 |
| 10 | turb vs B11 | 0.409447 |
Gráfica Turbidez en funcion de las distintas bandas
Código
import pandas as pd
import matplotlib.pyplot as plt
# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'turbidez']
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'turbidez']
# Cantidad de bandas
n = len(bandas)
# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3 # hasta 3 columnas por fila
columnas = 3
# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()
# Graficar cada banda
for i, banda in enumerate(bandas):
ax = axes[i]
ax.scatter(df[banda], df['turbidez'], alpha=0.6, s=10)
ax.set_title(f'Turbidez vs {banda}', fontsize=10)
ax.set_xlabel(banda)
ax.set_ylabel('Turbidez')
ax.grid(True)
# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
fig.delaxes(axes[j])
# Ajustar espacio
plt.tight_layout()
plt.show()Las mejores bandas para predecir turbidez (por correlación lineal) son:
B05 > B06 > B08 > B07 > B04
4.2 Prueba de correlación C2: logaritmo en la turbidez
Importamos nunpy para operar con funciones matemáticas
Creamos un nuevo DataFrame para aplicarle el logaritmo a la colomna turbidez.
Código
import numpy as np
Datos_turb_log = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
Datos_turb_log['turbidez'] = np.log(Datos_turb_log['turbidez'])
#Cambio el nombre la columna "turbidez" luego de aplicar el logaritmo
Datos_turb_log = Datos_turb_log.rename(columns={'turbidez': 'ln_turbidez'})
#Creo archivo csv para usarlo para entrenar modelos
Datos_turb_log.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv", index=False)
#print(Datos_turb_log.head())Calculamos r entre el ln_turb y cada banda, con la función .corr de pandas.
Código
bandas = [col for col in Datos_turb_log.columns if col.startswith('B')]
correlaciones = {}
for banda in bandas:
r = Datos_turb_log['ln_turbidez'].corr(Datos_turb_log[banda])
correlaciones[banda] = r
#Creamos un Data Frame.
df_correlaciones2 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones2 = df_correlaciones2.sort_values(by='r', ascending=False).reset_index(drop=True)
df_correlaciones2['Banda'] = "ln_turb vs " + df_correlaciones2['Banda'].astype(str)
df_correlaciones2.rename(columns={'Banda': 'Combinación 2'}, inplace=True)
df_correlaciones2.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv', index=False)
df_correlaciones2| Combinación 2 | r | |
|---|---|---|
| 0 | ln_turb vs B05 | 0.813367 |
| 1 | ln_turb vs B04 | 0.752422 |
| 2 | ln_turb vs B06 | 0.714446 |
| 3 | ln_turb vs B07 | 0.707941 |
| 4 | ln_turb vs B08 | 0.703117 |
| 5 | ln_turb vs B8A | 0.603662 |
| 6 | ln_turb vs B03 | 0.514314 |
| 7 | ln_turb vs B02 | 0.429268 |
| 8 | ln_turb vs B01 | 0.391705 |
| 9 | ln_turb vs B12 | 0.242050 |
| 10 | ln_turb vs B11 | 0.222510 |
Gráfica
Código
import pandas as pd
import matplotlib.pyplot as plt
# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv")
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']
# Cantidad de bandas
n = len(bandas)
# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3 # hasta 3 columnas por fila
columnas = 3
# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()
# Graficar cada banda
for i, banda in enumerate(bandas):
ax = axes[i]
ax.scatter(df[banda], df['ln_turbidez'], alpha=0.6, s=10)
ax.set_title(f'ln_turbidez vs {banda}', fontsize=10)
ax.set_xlabel(banda)
ax.set_ylabel('ln_turbidez')
ax.grid(True)
# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
fig.delaxes(axes[j])
# Ajustar espacio
plt.tight_layout()
plt.show()El orden general se mantiene, pero los valores de r bajan un poco al aplicar el logaritmo a la turbidez. Las mejores bandas siguen siendo: B05 > B04 > B06 > B07 > B08
Aunque r haya bajado ligeramente, esas bandas siguen siendo las más prometedoras predictoras lineales.
4.3 Prueba de correlación C3: logaritmo en las bandas
Anteriormente se estaba trabajando con una base de datos de reflectancias procesadas con Sen2Cor. Actualmente trabajamos con una nueva base de datos con las reflectancias corregidas por ACOLITE.
Acolite y Sen2Cor son dos herramientas de software utilizadas para la corrección atmosférica de imágenes satelitales. Acolite se destaca en ambientes acuáticos eutróficos, mientras que Sen2Cor es un procesador más general para la creación de productos Sentinel-2 de nivel 2A.
IMPORTANTE: Este cambio provocó errores en el código a partir de esta sección en adelante, el error se debe a que hay valores de reflectancia que son cero y negativos en la banda 12.
En las pruebas de correlación C3 y C4 aplicamos logaritmo a las bandas. Por defincion de logaritmo no ponemos tomar valores cero ni negativos. Para salvar este error, se procedió a filtrar filas con esos valores.
Creaamos un nuevo DataFrame para aplicarle el logaritmo a las bandas.
Código
import pandas as pd
# Leer los datos
Datos = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
# Guardar cantidad original de filas
filas_originales = Datos.shape[0]
# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
Datos_filtrado = Datos[(Datos > 0).all(axis=1)]
# Calcular cuántas filas se eliminaron
filas_finales = Datos_filtrado.shape[0]
filas_eliminadas = filas_originales - filas_finales
# Mostrar resultado
print(f"Se eliminaron {filas_eliminadas} filas que contenían ceros o valores negativos.")
# Si querés seguir trabajando con el DataFrame filtrado:
DatosC3_C4 = Datos_filtrado
#DatosC3_C4Se eliminaron 30 filas que contenían ceros o valores negativos.
Aplicamos logaritmo a las bandas
Código
import numpy as np
DatosC3 = DatosC3_C4.copy()
col = [col for col in DatosC3.columns if col.startswith('B')]
DatosC3[col] = np.log(DatosC3[col])
#Cambiamos en nombre las columnas, agremamos ln_ a cada columna
DatosC3.columns = ['ln_' + col for col in DatosC3.columns]
DatosC3.rename(columns={'ln_turbidez': 'turbidez'}, inplace=True)
#Creo archivo csv para usarlo para entrenar modelos
DatosC3.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv", index=False)
#DatosC3Calculamos r entre turb y ln de cada banda, con la función .corr de pandas.
Código
bandas_ln = [col for col in DatosC3.columns if col.startswith('ln_B')]
correlaciones = {}
for banda in bandas_ln:
r = DatosC3['turbidez'].corr(DatosC3[banda])
correlaciones[banda] = r
#Creamos un Data Frame.
df_correlaciones3 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones3 = df_correlaciones3.sort_values(by='r', ascending=False).reset_index(drop=True)
df_correlaciones3['Banda'] = "turb vs " + df_correlaciones3['Banda'].astype(str)
df_correlaciones3.rename(columns={'Banda': 'Combinación 3'}, inplace=True)
df_correlaciones3.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv', index=False)
df_correlaciones3| Combinación 3 | r | |
|---|---|---|
| 0 | turb vs ln_B05 | 0.734057 |
| 1 | turb vs ln_B06 | 0.717491 |
| 2 | turb vs ln_B07 | 0.706154 |
| 3 | turb vs ln_B04 | 0.699797 |
| 4 | turb vs ln_B08 | 0.698987 |
| 5 | turb vs ln_B8A | 0.641109 |
| 6 | turb vs ln_B03 | 0.598708 |
| 7 | turb vs ln_B02 | 0.543950 |
| 8 | turb vs ln_B01 | 0.533258 |
| 9 | turb vs ln_B11 | 0.352344 |
| 10 | turb vs ln_B12 | 0.351269 |
ln_B05 , ln_B06 , ln_B07 , ln_B08
Estas cuatro variables transformadas logarítmicamente se correlacionan de forma fuerte con la turbidez.
Gráfica
Código
import pandas as pd
import matplotlib.pyplot as plt
# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv")
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col.startswith('ln_B')]
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col.startswith('ln_B')]
# Cantidad de bandas
n = len(bandas)
# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3 # hasta 3 columnas por fila
columnas = 3
# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()
# Graficar cada banda
for i, banda in enumerate(bandas):
ax = axes[i]
ax.scatter(df[banda], df['turbidez'], alpha=0.6, s=10)
ax.set_title(f'Turbidez vs {banda}', fontsize=10)
ax.set_xlabel(banda)
ax.set_ylabel('Turbidez')
ax.grid(True)
# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
fig.delaxes(axes[j])
# Ajustar espacio
plt.tight_layout()
plt.show()4.4 Prueba de correlación C4: logaritmo en turbidez y en bandas
Primeramente verificamos que no hayan valores negativos y cero. Ya que por defincion de logaritno no los podremos utilizar.
Código
import pandas as pd
# Leer los datos
Datos = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
# Guardar cantidad original de filas
filas_originales = Datos.shape[0]
# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
Datos_filtrado = Datos[(Datos > 0).all(axis=1)]
# Calcular cuántas filas se eliminaron
filas_finales = Datos_filtrado.shape[0]
filas_eliminadas = filas_originales - filas_finales
# Mostrar resultado
print(f"Se eliminaron {filas_eliminadas} filas que contenían ceros o valores negativos.")
# Si querés seguir trabajando con el DataFrame filtrado:
DatosC3_C4 = Datos_filtrado
#DatosC3_C4Se eliminaron 30 filas que contenían ceros o valores negativos.
Creaamos un nuevo DataFrame para aplicarle el logaritmo a todas las columnas
Importamos nunpy para operar con funciones matemáticas
Código
import numpy as np
DatosC4 = np.log(DatosC3_C4)
#Cambiamos en nombre las columnas, agremamos ln_ a cada columna
DatosC4.columns = ['ln_' + col for col in DatosC4.columns]
#Creo archivo csv para usarlo para entrenar modelos
DatosC4.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv", index=False)
#DatosC4.head()Calculamos r entre el ln_turb y ln de cada banda, con la función .corr de pandas.
Código
bandas_ln = [col for col in DatosC4.columns if col.startswith('ln_B')]
correlaciones = {}
for banda in bandas_ln:
r = DatosC4['ln_turbidez'].corr(DatosC4[banda])
correlaciones[banda] = r
#Creamos un Data Frame.
df_correlaciones4 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones4 = df_correlaciones4.sort_values(by='r', ascending=False).reset_index(drop=True)
df_correlaciones4['Banda'] = "ln_turb vs " + df_correlaciones4['Banda'].astype(str)
df_correlaciones4.rename(columns={'Banda': 'Combinación 4'}, inplace=True)
df_correlaciones4.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv', index=False)
df_correlaciones4| Combinación 4 | r | |
|---|---|---|
| 0 | ln_turb vs ln_B05 | 0.844506 |
| 1 | ln_turb vs ln_B04 | 0.787566 |
| 2 | ln_turb vs ln_B06 | 0.770889 |
| 3 | ln_turb vs ln_B07 | 0.754235 |
| 4 | ln_turb vs ln_B08 | 0.740944 |
| 5 | ln_turb vs ln_B8A | 0.661797 |
| 6 | ln_turb vs ln_B03 | 0.533858 |
| 7 | ln_turb vs ln_B02 | 0.447893 |
| 8 | ln_turb vs ln_B01 | 0.397983 |
| 9 | ln_turb vs ln_B11 | 0.250618 |
| 10 | ln_turb vs ln_B12 | 0.223253 |
Gráfica
Código
import pandas as pd
import matplotlib.pyplot as plt
# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']
# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']
# Cantidad de bandas
n = len(bandas)
# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3 # hasta 3 columnas por fila
columnas = 3
# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()
# Graficar cada banda
for i, banda in enumerate(bandas):
ax = axes[i]
ax.scatter(df[banda], df['ln_turbidez'], alpha=0.6, s=10)
ax.set_title(f'ln_turbidez vs {banda}', fontsize=10)
ax.set_xlabel(banda)
ax.set_ylabel('ln_turbidez')
ax.grid(True)
# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
fig.delaxes(axes[j])
# Ajustar espacio
plt.tight_layout()
plt.show()B05, B04, B06, B07, B08 son los mejores candidatos como predictores en un modelo de regresión. Tienen una correlación fuerte con la turbidez, por encima de 0.78.
Comparación con combinaciones anteriores:
| Combinación | Variables transformadas | Mejor r |
|---|---|---|
| 1 | turbidez vs Banda | B05 (0.88) |
| 2 | ln(turb) vs Banda | B05 (0.80) |
| 3 | turbidez vs ln(Banda) | B05 (0.72) |
| 4 | ln(turb) vs ln(Banda) | B05 (0.87) |
5 Seleción de variables
Se resume en una tabla las pruebas de correlación realizadas.
Código
import pandas as pd
C1 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv")
C2= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv")
C3 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv")
C4 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv")
C_comb = pd.concat([C1, C2 ,C3, C4], axis=1)
C_comb| Combinación 1 | r | Combinación 2 | r | Combinación 3 | r | Combinación 4 | r | |
|---|---|---|---|---|---|---|---|---|
| 0 | turb vs B05 | 0.887299 | ln_turb vs B05 | 0.813367 | turb vs ln_B05 | 0.734057 | ln_turb vs ln_B05 | 0.844506 |
| 1 | turb vs B06 | 0.876827 | ln_turb vs B04 | 0.752422 | turb vs ln_B06 | 0.717491 | ln_turb vs ln_B04 | 0.787566 |
| 2 | turb vs B08 | 0.869149 | ln_turb vs B06 | 0.714446 | turb vs ln_B07 | 0.706154 | ln_turb vs ln_B06 | 0.770889 |
| 3 | turb vs B07 | 0.867741 | ln_turb vs B07 | 0.707941 | turb vs ln_B04 | 0.699797 | ln_turb vs ln_B07 | 0.754235 |
| 4 | turb vs B04 | 0.827019 | ln_turb vs B08 | 0.703117 | turb vs ln_B08 | 0.698987 | ln_turb vs ln_B08 | 0.740944 |
| 5 | turb vs B8A | 0.788088 | ln_turb vs B8A | 0.603662 | turb vs ln_B8A | 0.641109 | ln_turb vs ln_B8A | 0.661797 |
| 6 | turb vs B03 | 0.673358 | ln_turb vs B03 | 0.514314 | turb vs ln_B03 | 0.598708 | ln_turb vs ln_B03 | 0.533858 |
| 7 | turb vs B02 | 0.614315 | ln_turb vs B02 | 0.429268 | turb vs ln_B02 | 0.543950 | ln_turb vs ln_B02 | 0.447893 |
| 8 | turb vs B01 | 0.611305 | ln_turb vs B01 | 0.391705 | turb vs ln_B01 | 0.533258 | ln_turb vs ln_B01 | 0.397983 |
| 9 | turb vs B12 | 0.429200 | ln_turb vs B12 | 0.242050 | turb vs ln_B11 | 0.352344 | ln_turb vs ln_B11 | 0.250618 |
| 10 | turb vs B11 | 0.409447 | ln_turb vs B11 | 0.222510 | turb vs ln_B12 | 0.351269 | ln_turb vs ln_B12 | 0.223253 |
A partir del análisis con las correlaciones entre turbidez (y su logaritmo) y distintas bandas (y sus logaritmos). Se puede observar que en las combinaciones 1 y 4 las correlaciones son mas altas respecto a C2 y C3. Por que que se seguirá el análisis con estas combinaciones:
• Combinación 1: turb vs bandas
• Combinación 4: ln_turb vs ln_bandas
Se entrenarán modelos con ambas combinaciones por separado y se evalurá su desempeño.
Todos los modelos serán sometidos a bootstrapping con enfoque out-of-bag (OOB):
→ Donde se tomará el 75% de los datos con reemplazo para entrenamiento; → Luego se usa el 25% restante (no muestreados) como conjunto de prueba (test).
Proponemos avanzar hacia un modelo multivariable, utilizando el AIC [11] como criterio de seleccion de variables, este criterio es muy conveniente para encontrar el modelo de mejor balance entre ajuste y complejidad.
¿Qué hace el AIC?
→ Penaliza modelos con demasiadas variables para evitar sobreajuste.
→ Se busca minimizar el AIC: cuanto más bajo, mejor el modelo.
→ No requiere que las variables tengan alta correlación individual, pero sí que mejoren el modelo conjunto.
AIC=n⋅ln(RMSE²)+2k
5.1 Combinación 1: turb vs bandas
Código
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from utils import run_bootstrap_linear_regression_analysis # Asegurate de que esta función devuelve lo necesario
# Leer los datos
Datos = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
y = Datos['turbidez']
# Variables
variables_fijas = ['B05']
variables_a_agregar = ['B06', 'B08', 'B07', 'B04', 'B8A', 'B03', 'B02', 'B01', 'B12', 'B11']
# Resultados
resultados = []
# Configuración de bootstrapping
n_iteraciones_bootstrap = 1000
# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
variables_usadas = variables_fijas + variables_a_agregar[:i]
X = Datos[variables_usadas].values
# División entrenamiento/test
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, shuffle=True, random_state=42
)
# Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = run_bootstrap_linear_regression_analysis(
X_train, y_train.to_numpy(), n_iteraciones=n_iteraciones_bootstrap)
# Modelo final con coeficientes promedio
modelo_final = LinearRegression()
modelo_final.coef_ = coef_prom
modelo_final.intercept_ = intercept_prom
# Predicción sobre testeo
y_pred = modelo_final.predict(X_test)
r2 = r2_score(y_test, y_pred)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
# R² ajustado
n_obs = len(y_test)
n_vars = X_test.shape[1]
r2_ajustado = 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)
# AIC
residuals = y_test - y_pred
rss = np.sum(residuals ** 2)
k = X_test.shape[1] + 1 # +1 por el intercepto
aic = n_obs * np.log(rss / n_obs) + 2 * k
# Guardar resultados
resultados.append({
"variables": ", ".join(variables_usadas),
"num_variables": len(variables_usadas),
"R²_train (bootstrap)": r2_train_boot,
"RMSE_train (bootstrap)": rmse_train_boot,
"R²_test": r2,
"R²-ajustado": r2_ajustado,
"RMSE_test": rmse,
"AIC": aic
})
# Convertir a DataFrame
df_resultados = pd.DataFrame(resultados)
df_resultados| variables | num_variables | R²_train (bootstrap) | RMSE_train (bootstrap) | R²_test | R²-ajustado | RMSE_test | AIC | |
|---|---|---|---|---|---|---|---|---|
| 0 | B05 | 1 | 0.817414 | 35.380587 | 0.301995 | 0.258369 | 33.294120 | 130.193709 |
| 1 | B05, B06 | 2 | 0.818655 | 35.260133 | 0.333426 | 0.244549 | 32.535871 | 131.364356 |
| 2 | B05, B06, B08 | 3 | 0.825575 | 34.580912 | 0.278433 | 0.123812 | 33.851391 | 134.791284 |
| 3 | B05, B06, B08, B07 | 4 | 0.842854 | 32.823338 | 0.410136 | 0.228639 | 30.606533 | 133.163686 |
| 4 | B05, B06, B08, B07, B04 | 5 | 0.928892 | 22.079519 | 0.502045 | 0.294564 | 28.121155 | 132.114797 |
| 5 | B05, B06, B08, B07, B04, B8A | 6 | 0.943200 | 19.733626 | 0.762204 | 0.632497 | 19.433035 | 120.811081 |
| 6 | B05, B06, B08, B07, B04, B8A, B03 | 7 | 0.947473 | 18.976730 | 0.715553 | 0.516439 | 21.253927 | 126.035501 |
| 7 | B05, B06, B08, B07, B04, B8A, B03, B02 | 8 | 0.971103 | 14.075346 | 0.860893 | 0.737242 | 14.863235 | 115.160065 |
| 8 | B05, B06, B08, B07, B04, B8A, B03, B02, B01 | 9 | 0.976241 | 12.762692 | 0.902274 | 0.792332 | 12.457859 | 110.804661 |
| 9 | B05, B06, B08, B07, B04, B8A, B03, B02, B01, B12 | 10 | 0.979273 | 11.920504 | 0.925210 | 0.818367 | 10.898344 | 107.989991 |
| 10 | B05, B06, B08, B07, B04, B8A, B03, B02, B01, B... | 11 | 0.980290 | 11.624634 | 0.917927 | 0.767460 | 11.416642 | 111.662596 |
Gráfico de AIC
Código
aic_scores = [r["AIC"] for r in resultados]
n_vars = [r["num_variables"] for r in resultados]
plt.plot(n_vars, aic_scores, marker='o', color='purple')
plt.title('AIC vs Número de Variables')
plt.xlabel('Número de variables')
plt.ylabel('AIC')
plt.grid(True)
plt.show()¿Por qué el AIC baja al agregar más variables, si se supone que penaliza la complejidad?
AIC=n⋅ln(RMSE²)+2k
donde:
n = número de observaciones k = número de parámetros (incluye el intercepto) ln(RMSE²) representa el ajuste del modelo
AIC no busca modelos simples en sí, sino la mejor compensación entre ajuste y complejidad.
A pesar de agregar más variables, el error (RMSE) está disminuyendo mucho más rápido de lo que penaliza el AIC con +2k
5.2 Combinación 4: ln_turb vs ln_bandas
Código
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from utils import run_bootstrap_linear_regression_analysis # Asegurate de que esta función devuelve lo necesario
# Leer los datos
Datos = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
y = Datos['ln_turbidez']
# Variables
variables_fijas = ['ln_B05']
variables_a_agregar = ['ln_B04', 'ln_B06', 'ln_B07', 'ln_B08', 'ln_B8A', 'ln_B03', 'ln_B02', 'ln_B01','ln_B11' ,'ln_B12']
# Resultados
resultados = []
# Configuración de bootstrapping
n_iteraciones_bootstrap = 1000
# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
variables_usadas = variables_fijas + variables_a_agregar[:i]
X = Datos[variables_usadas].values
# División entrenamiento/test
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, shuffle=True, random_state=42
)
# Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = run_bootstrap_linear_regression_analysis(
X_train, y_train.to_numpy(), n_iteraciones=n_iteraciones_bootstrap)
# Modelo final con coeficientes promedio
modelo_final = LinearRegression()
modelo_final.coef_ = coef_prom
modelo_final.intercept_ = intercept_prom
# Predicción sobre testeo
y_pred = modelo_final.predict(X_test)
r2 = r2_score(y_test, y_pred)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
# Cálculo R² ajustado
n_obs = len(y_test)
n_vars = X_test.shape[1]
if n_obs - n_vars - 1 != 0:
r2_ajustado = 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)
else:
r2_ajustado = np.nan # o podrías usar 0.0 si preferís
# AIC
residuals = y_test - y_pred
rss = np.sum(residuals ** 2)
k = X_test.shape[1] + 1 # +1 por el intercepto
aic = n_obs * np.log(rss / n_obs) + 2 * k
# Guardar resultados
resultados.append({
"variables": ", ".join(variables_usadas),
"num_variables": len(variables_usadas),
"R²_train (bootstrap)": r2_train_boot,
"RMSE_train (bootstrap)": rmse_train_boot,
"R²_test": r2,
"R²-ajustado": r2_ajustado,
"RMSE_test": rmse,
"AIC": aic
})
# Convertir a DataFrame
df_resultados = pd.DataFrame(resultados)
df_resultados| variables | num_variables | R²_train (bootstrap) | RMSE_train (bootstrap) | R²_test | R²-ajustado | RMSE_test | AIC | |
|---|---|---|---|---|---|---|---|---|
| 0 | ln_B05 | 1 | 0.683149 | 0.731301 | 0.772278 | 0.746976 | 0.600810 | -7.208492 |
| 1 | ln_B05, ln_B04 | 2 | 0.814223 | 0.559970 | 0.764594 | 0.705742 | 0.610863 | -4.843421 |
| 2 | ln_B05, ln_B04, ln_B06 | 3 | 0.911953 | 0.385503 | 0.756236 | 0.651765 | 0.621613 | -2.459633 |
| 3 | ln_B05, ln_B04, ln_B06, ln_B07 | 4 | 0.909160 | 0.391569 | 0.732322 | 0.553871 | 0.651390 | 0.569765 |
| 4 | ln_B05, ln_B04, ln_B06, ln_B07, ln_B08 | 5 | 0.908863 | 0.392208 | 0.731631 | 0.463261 | 0.652231 | 2.598159 |
| 5 | ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A | 6 | 0.944083 | 0.307215 | 0.559205 | -0.101988 | 0.835898 | 10.056536 |
| 6 | ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... | 7 | 0.944316 | 0.306573 | 0.579665 | -0.401118 | 0.816268 | 11.533732 |
| 7 | ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... | 8 | 0.946340 | 0.300950 | 0.639149 | -0.804256 | 0.756309 | 11.855272 |
| 8 | ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... | 9 | 0.946290 | 0.301090 | 0.727218 | -1.727822 | 0.657571 | 10.777558 |
| 9 | ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... | 10 | 0.947751 | 0.296966 | 0.710570 | NaN | 0.677340 | 13.429195 |
| 10 | ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... | 11 | 0.951308 | 0.286679 | 0.696386 | 4.036139 | 0.693738 | 15.955473 |
Gráfico de AIC
Código
aic_scores = [r["AIC"] for r in resultados]
n_vars = [r["num_variables"] for r in resultados]
plt.plot(n_vars, aic_scores, marker='o', color='purple')
plt.title('AIC vs Número de Variables')
plt.xlabel('Número de variables')
plt.ylabel('AIC')
plt.grid(True)
plt.show()6 Entrenamiento de Modelos
6.1 Modelos sin logaritmo (turb vs bandas)
6.1.1 Modelo con B05
Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Variables predictoras
variables = ['B05']
X_completo = df[variables].values
y_completo = df['turbidez'].values
# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
X_completo, y_completo, test_size=0.25, random_state=42
)
# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)
# ------------------ Guardar coeficientes del modelo B05 ------------------ #
coef_B05 = coef_prom[0] # coeficiente de la banda B05
intercept_B05 = intercept_prom # intercepto
# Ecuación en LaTeX (solo para mostrar)
ecuacion_latexB05 = f"turbidez = {coef_B05:.4f} \\cdot B05 + {intercept_B05:.4f}"
from IPython.display import display, Math
display(Math(ecuacion_latexB05))
# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom
# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)
# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))
# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')
# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")
# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
# ------------------ Coeficientes ------------------ #
# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
'B05': 'B05'
}
# Construir términos con nombres renombrados
terminos_latex = [
f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
for var, coef in zip(variables, coef_prom)
]
# Construir la ecuación completa en LaTeX
ecuacion_latexB05 = " + ".join(terminos_latex)
ecuacion_latexB05 = f"turbidez = {ecuacion_latexB05} + {intercept_prom:.4f}"
display(Math(ecuacion_latexB05))\(\displaystyle turbidez = 1286.4244 \cdot B05 + -39.8511\)
\(\displaystyle turbidez = 1286.4244 \cdot B05 + -39.8511\)
6.1.2 Modelo con 6 bandas
Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Variables predictoras
variables = ['B05', 'B06', 'B08', 'B07', 'B04', 'B8A']
X_completo = df[variables].values
y_completo = df['turbidez'].values
# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
X_completo, y_completo, test_size=0.25, random_state=42
)
# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)
# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom
# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)
# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))
# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')
# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")
# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
# ------------------ Coeficientes ------------------ #
# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
'B05': 'B05',
'B04': 'B04',
'B06': 'B06',
'B07': 'B07',
'B8A': 'B8A',
'B08': 'B08'
}
# Construir términos con nombres renombrados
terminos_latex = [
f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
for var, coef in zip(variables, coef_prom)
]
# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"turbidez = {ecuacion_latex} + {intercept_prom:.4f}"
from IPython.display import display, Math
display(Math(ecuacion_latex))\(\displaystyle turbidez = 400.9979 \cdot B05 + 5763.2933 \cdot B06 + 1707.9858 \cdot B08 + -1647.7470 \cdot B07 + -1141.7309 \cdot B04 + -3926.9242 \cdot B8A + 2.4209\)
6.2 Modelos aplicando logaritmo en turbidez y en bandas
6.2.1 Modelo con 1 variable (ln_B05)
Escala logarítmica
Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
# Variables predictoras
variables = ['ln_B05']
X_completo = df[variables].values
y_completo = df['ln_turbidez'].values
# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
X_completo, y_completo, test_size=0.25, random_state=42
)
# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)
# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom
# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)
# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))
# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')
# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")
# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
# ------------------ Coeficientes ------------------ #
# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
'ln_B05': 'ln(B05)',
'ln_B04': 'ln(B04)',
'ln_B06': 'ln(B06)'
}
# Construir términos con nombres renombrados
terminos_latex = [
f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
for var, coef in zip(variables, coef_prom)
]
# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"
from IPython.display import display, Math
display(Math(ecuacion_latex))\(\displaystyle \ln(\text{turbidez}) = 1.5714 \cdot ln(B05) + 7.6043\)
Escala real
Código
# Transformo a escala real
y_train_pred_real = np.exp(y_train_pred)
y_test_pred_real = np.exp(y_test_pred)
y_train_real = np.exp(y_train)
y_test_real = np.exp(y_test)
rmse_test_real = np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
r2_test_real = r2_score(y_test_real, y_test_pred_real)
#Grafico
plt.figure(figsize=(4.5, 4.5))
# Entrenamiento (en escala real)
plt.scatter(y_train_real, y_train_pred_real, color="#9D50A6", alpha=0.5, label="Entrenamiento")
# Test (en escala real)
plt.scatter(y_test_real, y_test_pred_real, color="red", alpha=0.7, label="Test")
# Línea ideal
min_val = min(y_train_real.min(), y_test_real.min())
max_val = max(y_train_real.max(), y_test_real.max())
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")
# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Regresión lineal en escala real\n"
f"R² testeo: {r2_test_real:.4f}, RMSE: {rmse_test_real:.4f} NTU"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()6.2.2 Modelo con 3 variables
Escala Logarítmica
Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
# Variables predictoras
variables = ['ln_B05','ln_B04','ln_B06']
X_completo = df[variables].values
y_completo = df['ln_turbidez'].values
# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
X_completo, y_completo, test_size=0.25, random_state=42
)
# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)
# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom
# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)
# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))
# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')
# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")
# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
# ------------------ Coeficientes ------------------ #
# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
'ln_B05': 'ln(B05)',
'ln_B04': 'ln(B04)',
'ln_B06': 'ln(B06)'
}
# Construir términos con nombres renombrados
terminos_latex = [
f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
for var, coef in zip(variables, coef_prom)
]
# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"
from IPython.display import display, Math
display(Math(ecuacion_latex))
# Ecuacion real guardada como string
# Exponenciar el intercepto
intercepto_escala_real = np.exp(intercept_prom)
# Armar términos con potencias
terminos_real = [
f"{var[3:]}^{{{coef:.4f}}}" # Elimina "ln_" del nombre de la banda
for var, coef in zip(variables, coef_prom)
]
# Construir la ecuación en LaTeX como string reutilizable
ecuacion_real_latex = f"\\text{{turbidez}} = {intercepto_escala_real:.4f} \\cdot " + " \\cdot ".join(terminos_real)
# (Opcional) Mostrar
from IPython.display import Math, display
display(Math(ecuacion_real_latex))
# Guardar en otra variable si querés reutilizar luego como string sin formato LaTeX
ecuacion_turb = f"turbidez = {intercepto_escala_real:.4f} * " + " * ".join([
f"{var[3:]}^{coef:.4f}" for var, coef in zip(variables, coef_prom)
])\(\displaystyle \ln(\text{turbidez}) = 8.4385 \cdot ln(B05) + -5.4953 \cdot ln(B04) + -1.8635 \cdot ln(B06) + 6.1811\)
\(\displaystyle \text{turbidez} = 483.5190 \cdot B05^{8.4385} \cdot B04^{-5.4953} \cdot B06^{-1.8635}\)
Escala Real
Código
# Transformo a escala real
y_train_pred_real = np.exp(y_train_pred)
y_test_pred_real = np.exp(y_test_pred)
y_train_real = np.exp(y_train)
y_test_real = np.exp(y_test)
rmse_test_real = np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
r2_test_real = r2_score(y_test_real, y_test_pred_real)
#Grafico
plt.figure(figsize=(4.5, 4.5))
# Entrenamiento (en escala real)
plt.scatter(y_train_real, y_train_pred_real, color="#9D50A6", alpha=0.5, label="Entrenamiento")
# Test (en escala real)
plt.scatter(y_test_real, y_test_pred_real, color="red", alpha=0.7, label="Test")
# Línea ideal
min_val = min(y_train_real.min(), y_test_real.min())
max_val = max(y_train_real.max(), y_test_real.max())
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")
# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Regresión lineal en escala real\n"
f"R² testeo: {r2_test_real:.4f}, RMSE: {rmse_test_real:.4f} NTU"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()7 Procesamiento de productos satelitales
Para el procesamiento se utilizará la librería rasterio en python.
Tabla de referencia para seleccionar bandas en rasterio
| Nombre de banda en archivo .tif | Banda de Sentinel-2 | Número para la lectura en rasterio |
|---|---|---|
| Banda 1 | B01 | 1 |
| Banda 2 | B02 | 2 |
| Banda 3 | B03 | 3 |
| Banda 4 | B04 | 4 |
| Banda 5 | B05 | 5 |
| Banda 6 | B06 | 6 |
| Banda 7 | B07 | 7 |
| Banda 8 | B08 | 8 |
| Banda 9 | B8A | 9 |
| Banda 10 | B11 | 10 |
| Banda 11 | B12 | 11 |
Lectura de imágen raster y visualización del área de estudio, con el objetivo de verificar que el código está funcionando y de que se leen correctamente los archivos.
Código
import rasterio as rs
import matplotlib.pyplot as plt
import numpy as np
ruta = r"D:\GIT\Turbidez\recorte_acolite\2025-06-09.tif"
# Abrir el archivo
with rs.open(ruta) as raster:
#Definimos las bandas
#Bandas RGB para ver área de estudio
B02 = raster.read(2) #Azul
B03 = raster.read(3) #Verde
B04 = raster.read(4) #Rojo
#Para calcular NDWI con B03 y B11
B11 = raster.read(10)
#para el modelo
B05 = raster.read(5)
# Composión de imagen
rgb = np.stack([B04, B03, B02], axis=-1)
# Mejorar constraste
def stretch(banda):
p2, p98 = np.percentile(banda, (2, 98))
return np.clip((banda - p2) / (p98 - p2 + 1e-6), 0, 1)
rgb_stretched = np.stack([stretch(B04), stretch(B03), stretch(B02)], axis=-1)
# Gráfico
plt.figure(figsize=(4.5, 4.5))
plt.imshow(rgb_stretched)
plt.title("Río Paraná")
plt.axis("off")
plt.tight_layout()
plt.show()Se utilizrá el índice NDWI para recortar la imágen satelital, y poder trabajar únicamente con los píxeles de agua.
\(NDWI=(B03-B11)/(B03+B11)\)
Código
import numpy as np
import matplotlib.pyplot as plt
from skimage.filters import threshold_otsu
# NDWI
ndwi = (B03 - B11) / (B03 + B11 + 1e-6)
# Calcular umbral automático con Otsu
umbral = threshold_otsu(ndwi[np.isfinite(ndwi)]) # descartamos NaN e inf
# Binarización: agua = 1, no agua = 0
ndwi_bin = np.where(ndwi >= umbral, 1, 0)
# Gráfico máscara de agua
fig, ax = plt.subplots(figsize=(4.5, 4.5))
img = ax.imshow(ndwi_bin, cmap='Greys')
ax.set_title(f"NDWI - Otsu (umbral={umbral:.3f})")
cbar = fig.colorbar(img, ax=ax, fraction=0.046, pad=0.04)
cbar.set_label("Agua (1) / No Agua (0)")
ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
ax.axis("on")
plt.tight_layout()
plt.show()Luego del anális de las distantas pruebas de correlación con diferenctes bandas y combinaciones de ellas, optamos por tomar el modelo que utiliza la banda B05, para estimar la turbidez.
Aplicamos la ecuación hallada para ver la distribución de la turbidez en el río.
Código
import numpy as np
import matplotlib.pyplot as plt
import re
# --- Calcular turbidez directamente con variables guardadas ---
turbidez = coef_B05 * B05 + intercept_B05
# Forzar a que no existan valores negativos (los pone en 0)
turbidez = np.clip(turbidez, 0, None)
# --- Aplicar máscara NDWI (solo donde hay agua) ---
turbidez_agua = np.where(ndwi_bin == 1, turbidez, np.nan)
# --- Calcular percentiles para ajustar visualización ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])Mapa de turbidez
Código
import matplotlib.pyplot as plt
# Crear figura
fig, ax = plt.subplots(figsize=(4.5, 4.5))
# Mostrar imagen con coordenadas reales
img = ax.imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
# Título y colorbar
ax.set_title("Turbidez 2025-06-09")
cbar = fig.colorbar(img, ax=ax, fraction=0.046, pad=0.04)
cbar.set_label("Turbidez (NTU)")
ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
ax.axis("on")
plt.tight_layout()
plt.show()Para exportar tif
Código
import rasterio
import numpy as np
import os
import re
# Ruta de una banda original para copiarle la georreferenciación
ruta_base = r"D:\GIT\Turbidez\recorte_acolite\2025-06-09.tif"
# Carpeta donde guardar el resultado
carpeta_salida = r"D:\GIT\Turbidez\tif salidas"
# --- Extraer la fecha del nombre del archivo base ---
nombre_base = os.path.basename(ruta_base) # "2025-06-09.tif"
fecha = re.search(r"\d{4}-\d{2}-\d{2}", nombre_base).group() # "2025-06-09"
# Nombre final con fecha incluida
nombre_salida = f"turbidez_{fecha}.tif"
ruta_salida = os.path.join(carpeta_salida, nombre_salida)
# Abrimos la banda base para obtener el perfil (tamaño, proyección, etc.)
with rasterio.open(ruta_base) as src:
perfil = src.profile
# Ajustar perfil para nuestro resultado
perfil.update(
dtype=rasterio.float32,
count=1,
compress="lzw",
nodata=np.nan
)
# Guardar turbidez_agua en un GeoTIFF
with rasterio.open(ruta_salida, "w", **perfil) as dst:
dst.write(turbidez_agua.astype(rasterio.float32), 1)8 Lectura automática de tif y procesamiento para mapas de turbidez, cálculo de umbral de mácara de agua con método Otsu
Código
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
from skimage.filters import threshold_otsu
# Carpeta con archivos .tif
carpeta = r"D:\GIT\Turbidez\recorte_acolite"
# Buscar todos los .tif
tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"Se encontraron {len(tifs)} archivos .tif")
for ruta in tifs:
nombre = os.path.basename(ruta)
fecha = os.path.splitext(nombre)[0] # nombre sin extensión
print(f"\nProcesando: {nombre}")
with rs.open(ruta) as raster:
# Bandas necesarias
B03 = raster.read(3).astype(np.float32) # Verde
B05 = raster.read(5).astype(np.float32) # RedEdge (para ecuación turbidez)
B11 = raster.read(10).astype(np.float32) # SWIR
# --- NDWI ---
ndwi = (B03 - B11) / (B03 + B11 + 1e-6)
# --- Umbral Otsu ---
umbral = threshold_otsu(ndwi[np.isfinite(ndwi)])
print(f" Umbral NDWI (Otsu): {umbral:.3f}")
# --- Binarización: agua = 1, no agua = 0 ---
ndwi_bin = np.where(ndwi >= umbral, 1, 0)
# --- Turbidez ---
turbidez = coef_B05 * B05 + intercept_B05
turbidez = np.clip(turbidez, 0, None) # negativos a 0
turbidez_agua = np.where(ndwi_bin == 1, turbidez, np.nan)
# --- Calcular percentiles para ajustar visualización ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])
# --- Gráfico conjunto ---
fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))
# Panel 1: máscara de agua
im1 = axes[0].imshow(ndwi_bin, cmap="Greys")
axes[0].set_title(f"Máscara de agua - {fecha}")
axes[0].axis("off")
cbar1 = fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)
cbar1.set_label("Agua (1) / No Agua (0)")
# Panel 2: turbidez
im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
axes[1].set_title(f"Turbidez - {fecha}")
axes[1].axis("off")
cbar2 = fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04)
cbar2.set_label("Turbidez (NTU)")
plt.tight_layout()
plt.show()Se encontraron 7 archivos .tif
Procesando: 2023-05-11.tif
Umbral NDWI (Otsu): 0.308
Procesando: 2023-12-12.tif
Umbral NDWI (Otsu): -0.084
Procesando: 2024-05-20.tif
Umbral NDWI (Otsu): 37.771
Procesando: 2024-05-30.tif
Umbral NDWI (Otsu): -0.148
Procesando: 2024-12-16.tif
Umbral NDWI (Otsu): -138.005
Procesando: 2025-06-09.tif
Umbral NDWI (Otsu): 0.279
Procesando: 2025-09-12.tif
Umbral NDWI (Otsu): 1.153
Como el método Otsu no puede procesar correctamente algunas imágenes, se decidió utilizar el valor promedio del NDWI para calcular el umbral.
9 Conclusión para modelos entrenados con métodos de regresion lineal clásicos
Se seleccionará el mejor modelo lineal, este es el que utiliza la B05 de Sentinel-2 como variable predictora de la turbidez y será empleado para hacer la comparación con los métodos de aprendizaje automático próximamente. A continuación se resume la ecuación, métricas de desempeño y mapas de turbidez:
Modelo con B05 y sus métricas
Código
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Variables predictoras
variables = ['B05']
X_completo = df[variables].values
y_completo = df['turbidez'].values
# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
X_completo, y_completo, test_size=0.25, random_state=42
)
# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)
# ------------------ Guardar coeficientes del modelo B05 ------------------ #
coef_B05 = coef_prom[0] # coeficiente de la banda B05
intercept_B05 = intercept_prom # intercepto
# Ecuación en LaTeX (solo para mostrar)
ecuacion_latexB05 = f"turbidez = {coef_B05:.4f} \\cdot B05 + {intercept_B05:.4f}"
from IPython.display import display, Math
display(Math(ecuacion_latexB05))
# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom
# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)
# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))
# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')
# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')
# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")
# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1) # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
# ------------------ Coeficientes ------------------ #
# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
'B05': 'B05'
}
# Construir términos con nombres renombrados
terminos_latex = [
f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
for var, coef in zip(variables, coef_prom)
]
# Construir la ecuación completa en LaTeX
ecuacion_latexB05 = " + ".join(terminos_latex)
ecuacion_latexB05 = f"turbidez = {ecuacion_latexB05} + {intercept_prom:.4f}"
display(Math(ecuacion_latexB05))\(\displaystyle turbidez = 1294.7081 \cdot B05 + -40.1019\)
\(\displaystyle turbidez = 1294.7081 \cdot B05 + -40.1019\)
Mapas de turbidez
Código
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
from skimage.filters import threshold_otsu
# Carpeta con archivos .tif
carpeta = r"D:\GIT\Turbidez\recorte_acolite"
salida = os.path.join(carpeta, "turbidez_out") # subcarpeta para salidas
# Crear carpeta si no existe
os.makedirs(salida, exist_ok=True)
# Buscar todos los .tif
tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"Se encontraron {len(tifs)} archivos .tif")
for ruta in tifs:
nombre = os.path.basename(ruta)
fecha = os.path.splitext(nombre)[0] # nombre sin extensión
print(f"\nProcesando: {nombre}")
with rs.open(ruta) as raster:
# Bandas necesarias
B03 = raster.read(3).astype(np.float32) # Verde
B05 = raster.read(5).astype(np.float32) # RedEdge
B11 = raster.read(10).astype(np.float32) # SWIR
# Guardar perfil original
perfil = raster.profile
# --- NDWI ---
ndwi = (B03 - B11) / (B03 + B11 + 1e-6)
# Calcular el promedio (ignorando NaN si hay)
umbral = np.nanmean(ndwi)
# Binarizar con la media
mascara = (ndwi > umbral).astype(np.uint8)
# --- Turbidez ---
turbidez = coef_B05 * B05 + intercept_B05
turbidez = np.clip(turbidez, 0, None) # negativos a 0
turbidez_agua = np.where(mascara == 1, turbidez, np.nan)
# --- Exportar a GeoTIFF ---
salida_tif = os.path.join(salida, f"turbidez_{fecha}.tif")
perfil.update(
dtype=rs.float32,
count=1, # solo una banda (turbidez)
compress='lzw',
nodata=np.nan
)
with rs.open(salida_tif, 'w', **perfil) as dst:
dst.write(turbidez_agua.astype(np.float32), 1)
# --- Calcular percentiles para ajustar visualización ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])
# --- Gráfico conjunto ---
fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))
# Panel 1: máscara de agua
im1 = axes[0].imshow(mascara, cmap="Greys")
axes[0].set_title(f"Máscara de agua - {fecha}")
axes[0].axis("off")
cbar1 = fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)
cbar1.set_label("Agua (1) / No Agua (0)")
# Panel 2: turbidez
im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
axes[1].set_title(f"Turbidez - {fecha}")
axes[1].axis("off")
cbar2 = fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04)
cbar2.set_label("Turbidez (NTU)")
plt.tight_layout()
plt.show()Se encontraron 7 archivos .tif
Procesando: 2023-05-11.tif
Procesando: 2023-12-12.tif
Procesando: 2024-05-20.tif
Procesando: 2024-05-30.tif
Procesando: 2024-12-16.tif
Procesando: 2025-06-09.tif
Procesando: 2025-09-12.tif
10 Método de aprendizaje automático XGBoost
Código
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import pandas as pd
import numpy as np
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Definimos variables predictoras (bandas) y objetivo (turbidez medida)
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]
# Dividir en train/test
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, random_state=42)
# Modelo XGBoost
model = XGBRegressor(
n_estimators=500,
learning_rate=0.05,
max_depth=6,
subsample=0.8,
colsample_bytree=0.8,
random_state=42)
# Entrenar
model.fit(X_train, y_train)
# --- Predicciones ---
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# --- Métricas entrenamiento ---
r2_train = r2_score(y_train, y_train_pred)
rmse_train = np.sqrt(mean_squared_error(y_train, y_train_pred))
# --- Métricas testeo ---
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# --- Print resultados ---
print("=== Resultados Entrenamiento ===")
print(f"R² train: {r2_train:.3f}")
print(f"RMSE train: {rmse_train:.3f}")
print("\n=== Resultados Testeo ===")
print(f"R² test: {r2_test:.3f}")
print(f"RMSE test: {rmse_test:.3f}")
# Gráfico estimado vs real
import matplotlib.pyplot as plt
import numpy as np
plt.figure(figsize=(4.5, 4.5))
# --- Puntos entrenamiento ---
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5,
label="Datos de entrenamiento", marker='o')
# --- Puntos test ---
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7,
label="Datos de testeo", marker='^')
# --- Línea ideal 1:1 ---
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--',
color="#17A77E", lw=2, label="Línea ideal")
# --- Línea de tendencia sobre test ---
coef_linea = np.polyfit(y_test, y_test_pred, 1) # pendiente, intercepto
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black',
lw=2, label="Línea de tendencia (test)")
# --- Etiquetas ---
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Modelo XGBoost\n"
f"R² train: {r2_train:.3f}, RMSE: {rmse_train:.3f} | "
f"R² test: {r2_test:.3f}, RMSE: {rmse_test:.3f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()=== Resultados Entrenamiento ===
R² train: 1.000
RMSE train: 0.007
=== Resultados Testeo ===
R² test: 0.906
RMSE test: 12.208
10.1 Mapas con XGBoost
Código
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
from skimage.filters import threshold_otsu
import pandas as pd
# --- Configuración ---
carpeta = r"D:\GIT\Turbidez\recorte_acolite"
salida = os.path.join(carpeta, "turbidez_out_XBOOST")
os.makedirs(salida, exist_ok=True)
# Definir nombres de features (idénticos a entrenamiento)
feature_names = ["B01", "B02", "B03", "B04", "B05",
"B06", "B07", "B08", "B11", "B12", "B8A"]
# Buscar todos los .tif
tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"Se encontraron {len(tifs)} archivos .tif")
for ruta in tifs:
nombre = os.path.basename(ruta)
fecha = os.path.splitext(nombre)[0]
print(f"\nProcesando: {nombre}")
with rs.open(ruta) as raster:
# Leer bandas
B01 = raster.read(1).astype(np.float32)
B02 = raster.read(2).astype(np.float32)
B03 = raster.read(3).astype(np.float32)
B04 = raster.read(4).astype(np.float32)
B05 = raster.read(5).astype(np.float32)
B06 = raster.read(6).astype(np.float32)
B07 = raster.read(7).astype(np.float32)
B08 = raster.read(8).astype(np.float32)
B8A = raster.read(9).astype(np.float32)
B11 = raster.read(10).astype(np.float32)
B12 = raster.read(11).astype(np.float32)
perfil = raster.profile
# --- NDWI y máscara de agua ---
ndwi = (B03 - B11) / (B03 + B11 + 1e-6)
umbral = np.nanmean(ndwi)
mascara = (ndwi > umbral).astype(np.uint8)
# --- Crear array de features ---
filas, cols = B03.shape
X_img = np.stack([B01, B02, B03, B04, B05,
B06, B07, B08, B11, B12, B8A], axis=-1)
X_img_2d = X_img.reshape(-1, X_img.shape[-1])
# Convertir a DataFrame con nombres de columnas iguales al entrenamiento
X_img_2d_df = pd.DataFrame(X_img_2d, columns=feature_names)
# --- Predecir turbidez ---
turbidez_pred = model.predict(X_img_2d_df)
turbidez_pred = turbidez_pred.reshape(filas, cols)
# Aplicar máscara de agua
turbidez_agua = np.where(mascara == 1, turbidez_pred, np.nan)
# --- Exportar a GeoTIFF ---
salida_tif = os.path.join(salida, f"turbidez_{fecha}.tif")
perfil.update(dtype=rs.float32, count=1, compress='lzw', nodata=np.nan)
with rs.open(salida_tif, 'w', **perfil) as dst:
dst.write(turbidez_agua.astype(np.float32), 1)
# --- Percentiles para visualización ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])
# --- Gráfico ---
fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))
im1 = axes[0].imshow(mascara, cmap="Greys")
axes[0].set_title(f"Máscara de agua - {fecha}")
axes[0].axis("off")
cbar1 = fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)
cbar1.set_label("Agua (1) / No Agua (0)")
im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
axes[1].set_title(f"Turbidez (XGBoost) - {fecha}")
axes[1].axis("off")
cbar2 = fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04)
cbar2.set_label("Turbidez (NTU)")
plt.tight_layout()
plt.show()Se encontraron 7 archivos .tif
Procesando: 2023-05-11.tif
Procesando: 2023-12-12.tif
Procesando: 2024-05-20.tif
Procesando: 2024-05-30.tif
Procesando: 2024-12-16.tif
Procesando: 2025-06-09.tif
Procesando: 2025-09-12.tif
10.2 Hiperparámetros utilizados
• n_estimators=500 → número de arboles que se entrenan secuencialmente
• learning_rate=0.05 → indica que tanta aprende cada arbol
Es el “paso” con que el modelo corrige los errores. -Valores bajos (0.01–0.1) → aprendizaje más lento pero más estable. -Valores altos (>0.3) → puede sobreajustar rápido.
• max_depth=6 → Profundidad máxima de cada árbol.
Controla la complejidad de las interacciones entre las variables (bandas). Más profundo → aprende patrones complejos, pero también ruido.
• subsample=0.8 → Fracción de los datos usados para entrenar cada árbol.
Introduce aleatoriedad. Ayuda a reducir el sobreajuste, 0.8 significa que cada árbol usa el 80% de las filas del dataset.
• colsample_bytree=0.8 → Fracción de las variables (columnas) usadas para cada árbol. Igual que subsample, pero sobre las bandas.
Esto fuerza a los árboles a usar combinaciones distintas de bandas, lo que mejora la generalización.
0.8 = buena práctica para evitar que el modelo dependa demasiado de una sola banda.
Los hiperparámetros controlan el comportamiento del modelo, pero no se aprenden automáticamente del conjunto de datos. Por eso, necesitamos un método que busque las mejores combinaciones posibles para obtener el máximo rendimiento sin sobreajustar.
10.3 Optimización de hiperparámetros
Existen muchos métodos de optimizacion de hiperparámetros
- Grid Search (búsqueda en rejilla)
Prueba todas las combinaciones posibles de los valores que uno defina.
Por ejemplo, si se prueba 3 valores de max_depth y 3 de learning_rate, va a entrenar 3 × 3 = 9 modelos.
Usa validación cruzada (cross-validation) para evaluar cada combinación.
Ventajas
• Simple, exhaustivo y fácil de entender.
• Ideal cuando el número de combinaciones es pequeño.
Desventajas
• Muy lento si probás muchos parámetros o valores.
• No aprende de los resultados anteriores (prueba todo a ciegas).
- Random Search (búsqueda aleatoria)
En lugar de probar todas las combinaciones, elige al azar un número fijo de combinaciones.
Ejemplo: podrías probar 30 combinaciones elegidas aleatoriamente de un rango mucho más grande.
Ventajas
• Mucho más rápido que Grid Search.
• Puede explorar un espacio grande de parámetros.
• Buen punto medio entre precisión y costo computacional.
Desventajas
• Resultados algo aleatorios (aunque reproducibles con random_state).
• No garantiza probar la “mejor” combinación posible.
- Búsqueda bayesiana inteligente (Optuna / Hyperopt / scikit-optimize)
Usa inteligencia adaptativa: aprende de las pruebas anteriores y dirige las siguientes pruebas hacia las regiones más prometedoras.
Es una búsqueda inteligente, no aleatoria ni exhaustiva.
Ventajas
• Mucho más eficiente que Grid o Random Search.
• Encuentra buenos resultados con menos intentos.
• Ideal para modelos complejos como XGBoost.
Desventajas
• Requiere instalar librerías externas (como optuna).
• Más técnica, aunque muy potente.
10.4 Modelo XGBoost optimizado con Grid Search (Búsqueda por cuadrícula)
Código
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import pandas as pd
import numpy as np
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Definimos variables predictoras (bandas) y objetivo (turbidez medida)
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]
# Dividir en train/test
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, random_state=42)
# --- OPTIMIZACIÓN DE HIPERPARÁMETROS CON GRID SEARCH ---
from sklearn.model_selection import GridSearchCV
from xgboost import XGBRegressor
# Definimos las combinaciones a probar
param_grid = {
'max_depth': [4, 6, 8],
'learning_rate': [0.01, 0.05, 0.1],
'subsample': [0.7, 0.8, 1.0],
'colsample_bytree': [0.7, 0.8, 1.0],
'n_estimators': [300, 500, 700]
}
# Modelo base
xgb_model = XGBRegressor(random_state=42)
# Configuración del GridSearch
grid_search = GridSearchCV(
estimator=xgb_model,
param_grid=param_grid,
scoring='r2',
cv=5, # validación cruzada con 5 particiones
verbose=0,
n_jobs=1 )
# Entrenamos el GridSearch
grid_search.fit(X_train, y_train)
# Mostramos los mejores parámetros y resultados
print("\n=== RESULTADOS GRID SEARCH ===")
print("Mejores parámetros encontrados:")
import pandas as pd
from IPython.display import display
# Convertir los mejores parámetros en DataFrame
best_params_df = pd.DataFrame([grid_search.best_params_])
# Mostrar tabla
display(best_params_df)
# Guardamos el mejor modelo
model = grid_search.best_estimator_
# Entrenamos nuevamente con los mejores parámetros (opcional)
model.fit(X_train, y_train)
# --- Predicciones ---
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# --- Métricas entrenamiento ---
r2_train = r2_score(y_train, y_train_pred)
rmse_train = np.sqrt(mean_squared_error(y_train, y_train_pred))
# --- Métricas testeo ---
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# --- Print resultados ---
print("=== Resultados Entrenamiento ===")
print(f"R² train: {r2_train:.3f}")
print(f"RMSE train: {rmse_train:.3f}")
print("\n=== Resultados Testeo ===")
print(f"R² test: {r2_test:.3f}")
print(f"RMSE test: {rmse_test:.3f}")
# Gráfico estimado vs real
import matplotlib.pyplot as plt
import numpy as np
plt.figure(figsize=(4.5, 4.5))
# --- Puntos entrenamiento ---
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5,
label="Datos de entrenamiento", marker='o')
# --- Puntos test ---
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7,
label="Datos de testeo", marker='^')
# --- Línea ideal 1:1 ---
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--',
color="#17A77E", lw=2, label="Línea ideal")
# --- Línea de tendencia sobre test ---
coef_linea = np.polyfit(y_test, y_test_pred, 1) # pendiente, intercepto
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black',
lw=2, label="Línea de tendencia (test)")
# --- Etiquetas ---
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Modelo XGBoost (Optimizado con Grid Searh)\n"
f"R² train: {r2_train:.3f}, RMSE: {rmse_train:.3f} | "
f"R² test: {r2_test:.3f}, RMSE: {rmse_test:.3f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
=== RESULTADOS GRID SEARCH ===
Mejores parámetros encontrados:
| colsample_bytree | learning_rate | max_depth | n_estimators | subsample | |
|---|---|---|---|---|---|
| 0 | 0.7 | 0.05 | 4 | 300 | 0.7 |
=== Resultados Entrenamiento ===
R² train: 1.000
RMSE train: 0.595
=== Resultados Testeo ===
R² test: 0.904
RMSE test: 12.366
10.5 Modelo XGBoost optimizado con Optuna (Búsqueda bayesiana inteligente)
Código
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import pandas as pd
import numpy as np
# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Definimos variables predictoras (bandas) y objetivo (turbidez medida)
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]
# Dividir en train/test
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, random_state=42)
# --- OPTIMIZACIÓN DE HIPERPARÁMETROS CON OPTUNA ---
import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)
from xgboost import XGBRegressor
from sklearn.model_selection import cross_val_score
def objective(trial):
# Espacio de búsqueda de hiperparámetros
params = {
'max_depth': trial.suggest_int('max_depth', 3, 10),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3),
'subsample': trial.suggest_float('subsample', 0.6, 1.0),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
'n_estimators': trial.suggest_int('n_estimators', 200, 800),
'random_state': 42
}
# Modelo con los parámetros actuales
model = XGBRegressor(**params)
# Validación cruzada de 5 particiones
score = cross_val_score(model, X_train, y_train, scoring='r2', cv=5).mean()
return score
# Crear estudio de optimización
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=50, show_progress_bar=True)
# Mostrar resultados
print("\n=== RESULTADOS OPTUNA ===")
print("Mejores parámetros encontrados:")
import pandas as pd
from IPython.display import display
# Convertir los mejores parámetros a DataFrame
best_params_df = pd.DataFrame([study.best_params])
# Agregar la métrica objetivo (R²)
best_params_df["R2_promedio_CV"] = study.best_value
# Mostrar tabla
display(best_params_df)
# Entrenar el modelo final con los mejores parámetros
model = XGBRegressor(**study.best_params)
model.fit(X_train, y_train)
# --- Predicciones ---
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# --- Métricas entrenamiento ---
r2_train = r2_score(y_train, y_train_pred)
rmse_train = np.sqrt(mean_squared_error(y_train, y_train_pred))
# --- Métricas testeo ---
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))
# --- Print resultados ---
print("=== Resultados Entrenamiento ===")
print(f"R² train: {r2_train:.3f}")
print(f"RMSE train: {rmse_train:.3f}")
print("\n=== Resultados Testeo ===")
print(f"R² test: {r2_test:.3f}")
print(f"RMSE test: {rmse_test:.3f}")
# Gráfico estimado vs real
import matplotlib.pyplot as plt
import numpy as np
plt.figure(figsize=(4.5, 4.5))
# --- Puntos entrenamiento ---
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5,
label="Datos de entrenamiento", marker='o')
# --- Puntos test ---
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7,
label="Datos de testeo", marker='^')
# --- Línea ideal 1:1 ---
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--',
color="#17A77E", lw=2, label="Línea ideal")
# --- Línea de tendencia sobre test ---
coef_linea = np.polyfit(y_test, y_test_pred, 1) # pendiente, intercepto
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black',
lw=2, label="Línea de tendencia (test)")
# --- Etiquetas ---
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
f"Modelo XGBoost (Optimizado con Optuna)\n"
f"R² train: {r2_train:.3f}, RMSE: {rmse_train:.3f} | "
f"R² test: {r2_test:.3f}, RMSE: {rmse_test:.3f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()C:\Users\VICTOR\AppData\Local\Programs\Python\Python313\Lib\site-packages\tqdm\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
0%| | 0/50 [00:00<?, ?it/s]Best trial: 0. Best value: 0.897072: 0%| | 0/50 [00:00<?, ?it/s]Best trial: 0. Best value: 0.897072: 2%|▏ | 1/50 [00:00<00:17, 2.76it/s]Best trial: 0. Best value: 0.897072: 2%|▏ | 1/50 [00:00<00:17, 2.76it/s]Best trial: 0. Best value: 0.897072: 4%|▍ | 2/50 [00:00<00:14, 3.25it/s]Best trial: 2. Best value: 0.922563: 4%|▍ | 2/50 [00:01<00:14, 3.25it/s]Best trial: 2. Best value: 0.922563: 6%|▌ | 3/50 [00:01<00:16, 2.91it/s]Best trial: 2. Best value: 0.922563: 6%|▌ | 3/50 [00:01<00:16, 2.91it/s]Best trial: 2. Best value: 0.922563: 8%|▊ | 4/50 [00:01<00:14, 3.18it/s]Best trial: 2. Best value: 0.922563: 8%|▊ | 4/50 [00:01<00:14, 3.18it/s]Best trial: 2. Best value: 0.922563: 10%|█ | 5/50 [00:01<00:17, 2.53it/s]Best trial: 2. Best value: 0.922563: 10%|█ | 5/50 [00:02<00:17, 2.53it/s]Best trial: 2. Best value: 0.922563: 12%|█▏ | 6/50 [00:02<00:20, 2.12it/s]Best trial: 2. Best value: 0.922563: 12%|█▏ | 6/50 [00:03<00:20, 2.12it/s]Best trial: 2. Best value: 0.922563: 14%|█▍ | 7/50 [00:03<00:21, 1.96it/s]Best trial: 2. Best value: 0.922563: 14%|█▍ | 7/50 [00:03<00:21, 1.96it/s]Best trial: 2. Best value: 0.922563: 16%|█▌ | 8/50 [00:03<00:20, 2.02it/s]Best trial: 2. Best value: 0.922563: 16%|█▌ | 8/50 [00:03<00:20, 2.02it/s]Best trial: 2. Best value: 0.922563: 18%|█▊ | 9/50 [00:03<00:17, 2.36it/s]Best trial: 2. Best value: 0.922563: 18%|█▊ | 9/50 [00:04<00:17, 2.36it/s]Best trial: 2. Best value: 0.922563: 20%|██ | 10/50 [00:04<00:19, 2.08it/s]Best trial: 2. Best value: 0.922563: 20%|██ | 10/50 [00:05<00:19, 2.08it/s]Best trial: 2. Best value: 0.922563: 22%|██▏ | 11/50 [00:05<00:20, 1.88it/s]Best trial: 2. Best value: 0.922563: 22%|██▏ | 11/50 [00:05<00:20, 1.88it/s]Best trial: 2. Best value: 0.922563: 24%|██▍ | 12/50 [00:05<00:21, 1.75it/s]Best trial: 2. Best value: 0.922563: 24%|██▍ | 12/50 [00:06<00:21, 1.75it/s]Best trial: 2. Best value: 0.922563: 26%|██▌ | 13/50 [00:06<00:22, 1.64it/s]Best trial: 2. Best value: 0.922563: 26%|██▌ | 13/50 [00:07<00:22, 1.64it/s]Best trial: 2. Best value: 0.922563: 28%|██▊ | 14/50 [00:07<00:22, 1.59it/s]Best trial: 2. Best value: 0.922563: 28%|██▊ | 14/50 [00:07<00:22, 1.59it/s]Best trial: 2. Best value: 0.922563: 30%|███ | 15/50 [00:07<00:23, 1.48it/s]Best trial: 2. Best value: 0.922563: 30%|███ | 15/50 [00:08<00:23, 1.48it/s]Best trial: 2. Best value: 0.922563: 32%|███▏ | 16/50 [00:08<00:21, 1.57it/s]Best trial: 16. Best value: 0.923312: 32%|███▏ | 16/50 [00:08<00:21, 1.57it/s]Best trial: 16. Best value: 0.923312: 34%|███▍ | 17/50 [00:08<00:19, 1.67it/s]Best trial: 17. Best value: 0.933624: 34%|███▍ | 17/50 [00:09<00:19, 1.67it/s]Best trial: 17. Best value: 0.933624: 36%|███▌ | 18/50 [00:09<00:17, 1.82it/s]Best trial: 17. Best value: 0.933624: 36%|███▌ | 18/50 [00:09<00:17, 1.82it/s]Best trial: 17. Best value: 0.933624: 38%|███▊ | 19/50 [00:09<00:15, 2.03it/s]Best trial: 17. Best value: 0.933624: 38%|███▊ | 19/50 [00:10<00:15, 2.03it/s]Best trial: 17. Best value: 0.933624: 40%|████ | 20/50 [00:10<00:15, 1.96it/s]Best trial: 17. Best value: 0.933624: 40%|████ | 20/50 [00:10<00:15, 1.96it/s]Best trial: 17. Best value: 0.933624: 42%|████▏ | 21/50 [00:10<00:13, 2.11it/s]Best trial: 17. Best value: 0.933624: 42%|████▏ | 21/50 [00:10<00:13, 2.11it/s]Best trial: 17. Best value: 0.933624: 44%|████▍ | 22/50 [00:10<00:11, 2.39it/s]Best trial: 17. Best value: 0.933624: 44%|████▍ | 22/50 [00:11<00:11, 2.39it/s]Best trial: 17. Best value: 0.933624: 46%|████▌ | 23/50 [00:11<00:10, 2.46it/s]Best trial: 17. Best value: 0.933624: 46%|████▌ | 23/50 [00:11<00:10, 2.46it/s]Best trial: 17. Best value: 0.933624: 48%|████▊ | 24/50 [00:11<00:11, 2.34it/s]Best trial: 17. Best value: 0.933624: 48%|████▊ | 24/50 [00:12<00:11, 2.34it/s]Best trial: 17. Best value: 0.933624: 50%|█████ | 25/50 [00:12<00:10, 2.36it/s]Best trial: 17. Best value: 0.933624: 50%|█████ | 25/50 [00:12<00:10, 2.36it/s]Best trial: 17. Best value: 0.933624: 52%|█████▏ | 26/50 [00:12<00:11, 2.09it/s]Best trial: 17. Best value: 0.933624: 52%|█████▏ | 26/50 [00:13<00:11, 2.09it/s]Best trial: 17. Best value: 0.933624: 54%|█████▍ | 27/50 [00:13<00:11, 1.92it/s]Best trial: 17. Best value: 0.933624: 54%|█████▍ | 27/50 [00:13<00:11, 1.92it/s]Best trial: 17. Best value: 0.933624: 56%|█████▌ | 28/50 [00:13<00:11, 1.86it/s]Best trial: 17. Best value: 0.933624: 56%|█████▌ | 28/50 [00:14<00:11, 1.86it/s]Best trial: 17. Best value: 0.933624: 58%|█████▊ | 29/50 [00:14<00:11, 1.91it/s]Best trial: 17. Best value: 0.933624: 58%|█████▊ | 29/50 [00:14<00:11, 1.91it/s]Best trial: 17. Best value: 0.933624: 60%|██████ | 30/50 [00:14<00:10, 2.00it/s]Best trial: 17. Best value: 0.933624: 60%|██████ | 30/50 [00:15<00:10, 2.00it/s]Best trial: 17. Best value: 0.933624: 62%|██████▏ | 31/50 [00:15<00:10, 1.85it/s]Best trial: 17. Best value: 0.933624: 62%|██████▏ | 31/50 [00:16<00:10, 1.85it/s]Best trial: 17. Best value: 0.933624: 64%|██████▍ | 32/50 [00:16<00:09, 1.92it/s]Best trial: 17. Best value: 0.933624: 64%|██████▍ | 32/50 [00:16<00:09, 1.92it/s]Best trial: 17. Best value: 0.933624: 66%|██████▌ | 33/50 [00:16<00:09, 1.87it/s]Best trial: 17. Best value: 0.933624: 66%|██████▌ | 33/50 [00:17<00:09, 1.87it/s]Best trial: 17. Best value: 0.933624: 68%|██████▊ | 34/50 [00:17<00:08, 1.90it/s]Best trial: 17. Best value: 0.933624: 68%|██████▊ | 34/50 [00:17<00:08, 1.90it/s]Best trial: 17. Best value: 0.933624: 70%|███████ | 35/50 [00:17<00:07, 1.93it/s]Best trial: 17. Best value: 0.933624: 70%|███████ | 35/50 [00:18<00:07, 1.93it/s]Best trial: 17. Best value: 0.933624: 72%|███████▏ | 36/50 [00:18<00:07, 1.81it/s]Best trial: 17. Best value: 0.933624: 72%|███████▏ | 36/50 [00:18<00:07, 1.81it/s]Best trial: 17. Best value: 0.933624: 74%|███████▍ | 37/50 [00:18<00:06, 1.89it/s]Best trial: 17. Best value: 0.933624: 74%|███████▍ | 37/50 [00:18<00:06, 1.89it/s]Best trial: 17. Best value: 0.933624: 76%|███████▌ | 38/50 [00:18<00:05, 2.22it/s]Best trial: 17. Best value: 0.933624: 76%|███████▌ | 38/50 [00:19<00:05, 2.22it/s]Best trial: 17. Best value: 0.933624: 78%|███████▊ | 39/50 [00:19<00:05, 2.16it/s]Best trial: 17. Best value: 0.933624: 78%|███████▊ | 39/50 [00:19<00:05, 2.16it/s]Best trial: 17. Best value: 0.933624: 80%|████████ | 40/50 [00:19<00:04, 2.35it/s]Best trial: 17. Best value: 0.933624: 80%|████████ | 40/50 [00:20<00:04, 2.35it/s]Best trial: 17. Best value: 0.933624: 82%|████████▏ | 41/50 [00:20<00:03, 2.38it/s]Best trial: 17. Best value: 0.933624: 82%|████████▏ | 41/50 [00:20<00:03, 2.38it/s]Best trial: 17. Best value: 0.933624: 84%|████████▍ | 42/50 [00:20<00:03, 2.28it/s]Best trial: 17. Best value: 0.933624: 84%|████████▍ | 42/50 [00:21<00:03, 2.28it/s]Best trial: 17. Best value: 0.933624: 86%|████████▌ | 43/50 [00:21<00:03, 1.79it/s]Best trial: 17. Best value: 0.933624: 86%|████████▌ | 43/50 [00:21<00:03, 1.79it/s]Best trial: 17. Best value: 0.933624: 88%|████████▊ | 44/50 [00:21<00:03, 1.93it/s]Best trial: 17. Best value: 0.933624: 88%|████████▊ | 44/50 [00:22<00:03, 1.93it/s]Best trial: 17. Best value: 0.933624: 90%|█████████ | 45/50 [00:22<00:02, 1.97it/s]Best trial: 17. Best value: 0.933624: 90%|█████████ | 45/50 [00:23<00:02, 1.97it/s]Best trial: 17. Best value: 0.933624: 92%|█████████▏| 46/50 [00:23<00:02, 1.89it/s]Best trial: 17. Best value: 0.933624: 92%|█████████▏| 46/50 [00:23<00:02, 1.89it/s]Best trial: 17. Best value: 0.933624: 94%|█████████▍| 47/50 [00:23<00:01, 1.91it/s]Best trial: 17. Best value: 0.933624: 94%|█████████▍| 47/50 [00:24<00:01, 1.91it/s]Best trial: 17. Best value: 0.933624: 96%|█████████▌| 48/50 [00:24<00:01, 1.93it/s]Best trial: 17. Best value: 0.933624: 96%|█████████▌| 48/50 [00:24<00:01, 1.93it/s]Best trial: 17. Best value: 0.933624: 98%|█████████▊| 49/50 [00:24<00:00, 1.86it/s]Best trial: 17. Best value: 0.933624: 98%|█████████▊| 49/50 [00:25<00:00, 1.86it/s]Best trial: 17. Best value: 0.933624: 100%|██████████| 50/50 [00:25<00:00, 1.84it/s]Best trial: 17. Best value: 0.933624: 100%|██████████| 50/50 [00:25<00:00, 1.98it/s]
=== RESULTADOS OPTUNA ===
Mejores parámetros encontrados:
| max_depth | learning_rate | subsample | colsample_bytree | n_estimators | R2_promedio_CV | |
|---|---|---|---|---|---|---|
| 0 | 6 | 0.071912 | 0.672075 | 0.942038 | 318 | 0.933624 |
=== Resultados Entrenamiento ===
R² train: 1.000
RMSE train: 0.082
=== Resultados Testeo ===
R² test: 0.942
RMSE test: 9.568
10.6 Relevancia de Bandas espectrales en el modelo XGBOOST
La importancia de las bandas indica cuánto contribuye cada banda espectral (o combinación de ellas) al rendimiento del modelo. Esto permite saber cuáles longitudes de onda son más informativas para estimar la variable objetivo.
Interpretación:
Cuanto mayor sea la importancia, más influye esa banda en las divisiones del modelo.
No necesariamente significa correlación directa (puede ser no lineal).
Si una banda tiene muy baja importancia (≈0), puede eliminarse en futuras versiones del modelo para simplificarlo.
Cuando entrenás un modelo como XGBoost, el algoritmo aprende qué variables (bandas espectrales) son más útiles para predecir la turbidez.
• La importancia de las bandas indica cuánto aportó cada banda al modelo para mejorar sus predicciones.
En XGBoost, el atributo model.feature_importances_ mide la importancia promedio de cada variable dentro de todos los árboles del modelo.
XGBoost da la importancia relativa (no absoluta ni porcentual) basada en “weight”. Los valores se normalizan para que la suma de todos sea 1.0.
Tipo de importancia (importance_type) |
Qué mide exactamente | Interpretación práctica |
|---|---|---|
"weight" (por defecto) |
Número de veces que la variable fue usada para dividir nodos | Bandas usadas más veces para tomar decisiones |
"gain" |
Promedio de la mejora en la predicción al usar esa banda | Cuánto “valor” aporta cada vez que se usa |
"cover" |
Promedio de la cantidad de datos afectados por esas divisiones | Qué tan general es su influencia |
"total_gain" o "total_cover" |
Suma total de las métricas anteriores | Importancia acumulada |
Bandas de Sentinel-2
Importancia de cada banda dentro del modelo XGBoost
Código
# === IMPORTANCIA DE BANDAS (Feature Importances) ===
import pandas as pd
import matplotlib.pyplot as plt
# Obtener importancias del modelo entrenado
importancias = model.feature_importances_
# Crear DataFrame ordenado
importancia_df = pd.DataFrame({
'Banda': X.columns,
'Importancia': importancias
}).sort_values(by='Importancia', ascending=False)
display(importancia_df)
# Gráfico de barras ordenado
plt.figure(figsize=(8, 4))
plt.bar(importancia_df['Banda'], importancia_df['Importancia'], color='#2E86AB')
plt.title('Importancia de las bandas en el modelo XGBoost')
plt.xlabel('Banda')
plt.ylabel('Importancia relativa')
plt.xticks(rotation=45)
plt.grid(axis='y', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()| Banda | Importancia | |
|---|---|---|
| 4 | B05 | 0.719612 |
| 5 | B06 | 0.157507 |
| 0 | B01 | 0.052247 |
| 7 | B08 | 0.027199 |
| 1 | B02 | 0.017780 |
| 10 | B8A | 0.009735 |
| 8 | B11 | 0.005484 |
| 6 | B07 | 0.002954 |
| 3 | B04 | 0.002905 |
| 2 | B03 | 0.002719 |
| 9 | B12 | 0.001858 |
B05 explica la mayor parte del poder predictivo total del modelo.
Las bandas con valores cercanos a cero aportan muy poca información o son redundantes con otras.
Las bandas Red Edge (B05–B06-B07) se ubican justo en la transición entre el rojo y el infrarrojo cercano, donde la reflectancia del agua cambia mucho con la concentración de sedimentos. → Son excelentes indicadores de turbidez o material particulado suspendido.
Las bandas azules (B01–B02) todavía aportan algo, porque la turbidez afecta la dispersión del azul. → En aguas claras, el azul predomina; en aguas turbias, se atenúa.
Las bandas SWIR (B11, B12) y el rojo visible (B04) no aportan casi nada, ya que esas longitudes de onda no penetran el agua y son más útiles para suelo o vegetación.
10.7 Mapas XGBOOST con hiperparámetros optimizados
Código
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
import pandas as pd
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
import optuna
from skimage.filters import threshold_otsu
# --- CONFIGURACIÓN ---
carpeta = r"D:\GIT\Turbidez\recorte_acolite"
salida = os.path.join(carpeta, "turbidez_out_XGBOOST_optuna")
os.makedirs(salida, exist_ok=True)
# --- ENTRENAMIENTO DEL MODELO OPTIMIZADO ---
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
# Variables predictoras y objetivo
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)
optuna.logging.set_verbosity(optuna.logging.WARNING)
def objective(trial):
params = {
'max_depth': trial.suggest_int('max_depth', 3, 10),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3),
'subsample': trial.suggest_float('subsample', 0.6, 1.0),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
'n_estimators': trial.suggest_int('n_estimators', 200, 800),
'random_state': 42
}
model = XGBRegressor(**params)
score = cross_val_score(model, X_train, y_train, scoring='r2', cv=5).mean()
return score
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=50, show_progress_bar=False)
best_params = study.best_params
print("=== MEJORES PARÁMETROS ===")
print(best_params)
print(f"R² promedio CV: {study.best_value:.3f}")
# Entrenar modelo final
model = XGBRegressor(**best_params)
model.fit(X_train, y_train)
# Evaluar desempeño
y_test_pred = model.predict(X_test)
print(f"R² test: {r2_score(y_test, y_test_pred):.3f}")
print(f"RMSE test: {np.sqrt(mean_squared_error(y_test, y_test_pred)):.3f}")
# --- PROCESAR RASTERS ---
feature_names = ["B01", "B02", "B03", "B04", "B05",
"B06", "B07", "B08", "B11", "B12", "B8A"]
tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"\nSe encontraron {len(tifs)} archivos .tif")
for ruta in tifs:
nombre = os.path.basename(ruta)
fecha = os.path.splitext(nombre)[0]
print(f"\nProcesando: {nombre}")
with rs.open(ruta) as raster:
bandas = [raster.read(i).astype(np.float32) for i in range(1, 12)]
perfil = raster.profile
# --- NDWI y máscara de agua ---
B03, B11 = bandas[2], bandas[9]
ndwi = (B03 - B11) / (B03 + B11 + 1e-6)
umbral = np.nanmean(ndwi)
mascara = (ndwi > umbral).astype(np.uint8)
# --- Crear matriz de predicción ---
filas, cols = B03.shape
X_img = np.stack(bandas, axis=-1)
X_img_2d = X_img.reshape(-1, X_img.shape[-1])
X_img_df = pd.DataFrame(X_img_2d, columns=feature_names)
# --- Predecir turbidez ---
turbidez_pred = model.predict(X_img_df)
turbidez_pred = turbidez_pred.reshape(filas, cols)
# Aplicar máscara
turbidez_agua = np.where(mascara == 1, turbidez_pred, np.nan)
# --- Exportar GeoTIFF ---
salida_tif = os.path.join(salida, f"turbidez_{fecha}.tif")
perfil.update(dtype=rs.float32, count=1, compress='lzw', nodata=np.nan)
with rs.open(salida_tif, 'w', **perfil) as dst:
dst.write(turbidez_agua.astype(np.float32), 1)
# --- Visualización ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])
fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))
im1 = axes[0].imshow(mascara, cmap="Greys")
axes[0].set_title(f"Máscara de agua - {fecha}")
axes[0].axis("off")
fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04, label="Agua / No Agua")
im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
axes[1].set_title(f"Turbidez (XGBoost Optuna) - {fecha}")
axes[1].axis("off")
fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04, label="Turbidez (NTU)")
plt.tight_layout()
plt.show()=== MEJORES PARÁMETROS ===
{'max_depth': 3, 'learning_rate': 0.022172460819822562, 'subsample': 0.7605112208446342, 'colsample_bytree': 0.7651599773171959, 'n_estimators': 372}
R² promedio CV: 0.935
R² test: 0.902
RMSE test: 12.482
Se encontraron 7 archivos .tif
Procesando: 2023-05-11.tif
Procesando: 2023-12-12.tif
Procesando: 2024-05-20.tif
Procesando: 2024-05-30.tif
Procesando: 2024-12-16.tif
Procesando: 2025-06-09.tif
Procesando: 2025-09-12.tif